Sesión 05

Diseño de Experimentos Factoriales

Christian Vásquez-Velasco, Bach., M.Sc.(c)

InkaStats Academy

2023

[Curso: Métodos Estadísticos para Experimentos Agrícolas en R]

Instalar paquetes necesarios

# install.packages("devtools")
# devtools::install_github("emitanaka/edibble", force = T)
# devtools::install_github("emitanaka/deggust", force = T)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(readxl, agricolae, agricolaeplotr, car, tidyverse, PMCMRplus, outliers, nortest, mvtnorm, lmtest, ExpDes, edibble, gt,
               gtsummary, devtools, deggust, xlsx, desplot,
               ggResidpanel, fastGraph, gvlma, multcomp,
               phia, dvmisc)

Diseño completamente al azar en arreglo factorial con dos vías

Planeamiento

Crear un libro de campo con el paquete agricolae

trt <- c(3,4)

r <- 10

salida <- agricolae::design.ab(trt = trt,
                               r = r,
                               design = "crd",
                               serie = 3,
                               seed = 123,
                               kinds = "Super-Duper",
                               randomization = TRUE)
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Guardar el libro generado

write.table(salida$book,
            "books/crd2w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida$book,
           "books/crd2w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
agricolaeplotr::plot_design.factorial_crd(
  design = salida,
  factor_name = "A",
  ncols = 12,
  nrows = 10,
  reverse_y = TRUE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas")
agricolaeplotr::plot_design.factorial_crd(
  design = salida,
  factor_name = "B",
  ncols = 12,
  nrows = 10,
  reverse_y = TRUE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Columnas",
       y = "Filas")

Crear un libro de campo con el paquete edibble

menu_factorial()
design("Factorial Design") %>%
  set_units(unit = 126) %>%
  set_trts(trt1 = 7,
           trt2 = 2) %>%
  allot_trts(~unit) %>%
  assign_trts("random", seed = 383) %>%
  serve_table()
crd <- takeout(menu_factorial(trt = c(3,4),
                              r = 10,
                              seed = 123))
crd %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
crd2 <- design("Completely Randomised Design") %>%
  set_units(unit = 120) %>%
  set_trts(trt1 = c("Amazon", "Coari x Lame", "Coari x Yangambí"),
           trt2 = c(0, 30, 60, 90)) %>%
  allot_trts(trt1 + trt2 ~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
crd2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
deggust::autoplot(crd)
plot(crd)

Análisis de DCA en arreglo factorial con dos vías sin interacción entre factores

Importación de datos

archivos <- list.files(pattern = "datos dca 2w.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor)

Creación del modelo lineal

modelo.dca <- lm(Rendimiento ~ Variedad + Densidad, data = data)

Definición del modelo

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2\]

summary(modelo.dca)

Call:
lm(formula = Rendimiento ~ Variedad + Densidad, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
 -3.25  -1.75   0.00   1.25   3.25 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.2500     0.7703  10.710 5.61e-09 ***
Variedada2    2.5000     0.8895   2.811  0.01203 *  
Densidadb2   -3.5000     0.8895  -3.935  0.00107 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.989 on 17 degrees of freedom
Multiple R-squared:  0.579, Adjusted R-squared:  0.5295 
F-statistic: 11.69 on 2 and 17 DF,  p-value: 0.0006399

Verificación visual de los supuestos del modelo

performance::check_model(modelo.dca)
ggResidpanel::resid_panel(modelo.dca)
influence.measures(modelo.dca)
Influence measures of
     lm(formula = Rendimiento ~ Variedad + Densidad, data = data) :

    dfb.1_ dfb.Vrd2 dfb.Dns2   dffit cov.r  cook.d  hat inf
1  -0.8000   0.4619   0.4619 -0.8000 0.765 0.18478 0.15    
2   0.3998  -0.2308  -0.2308  0.3998 1.196 0.05358 0.15    
3  -0.0556   0.0321   0.0321 -0.0556 1.407 0.00109 0.15    
4  -0.0556   0.0321   0.0321 -0.0556 1.407 0.00109 0.15    
5   0.1675  -0.0967  -0.0967  0.1675 1.370 0.00984 0.15    
6   0.1746  -0.3024   0.3024  0.5238 1.068 0.08856 0.15    
7  -0.1333   0.2308  -0.2308 -0.3998 1.196 0.05358 0.15    
8   0.0939  -0.1626   0.1626  0.2817 1.299 0.02733 0.15    
9   0.0185  -0.0321   0.0321  0.0556 1.407 0.00109 0.15    
10 -0.0558   0.0967  -0.0967 -0.1675 1.370 0.00984 0.15    
11 -0.2187  -0.3788   0.3788 -0.6561 0.922 0.13230 0.15    
12  0.2667   0.4619  -0.4619  0.8000 0.765 0.18478 0.15    
13  0.0939   0.1626  -0.1626  0.2817 1.299 0.02733 0.15    
14 -0.1333  -0.2308   0.2308 -0.3998 1.196 0.05358 0.15    
15  0.0939   0.1626  -0.1626  0.2817 1.299 0.02733 0.15    
16 -0.2187   0.3788   0.3788  0.6561 0.922 0.13230 0.15    
17  0.1746  -0.3024  -0.3024 -0.5238 1.068 0.08856 0.15    
18  0.1746  -0.3024  -0.3024 -0.5238 1.068 0.08856 0.15    
19  0.0185  -0.0321  -0.0321 -0.0556 1.407 0.00109 0.15    
20 -0.0558   0.0967   0.0967  0.1675 1.370 0.00984 0.15    
influenceIndexPlot(modelo.dca)

Cumplimiento de supuestos del modelo lineal general

Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dca,
                 reps = 5000,
                 max.lag = 5)
 lag Autocorrelation D-W Statistic p-value
   1      -0.2239777      2.282528  0.8296
   2      -0.3847584      2.557621  0.2228
   3       0.1830855      1.345725  0.2452
   4       0.2825279      1.070632  0.1568
   5      -0.5845725      2.684015  0.0040
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dca, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dca
DW = 2.2825, p-value = 0.8366
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dca))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dca)
W = 0.97395, p-value = 0.835
ad.test(rstudent(modelo.dca))

    Anderson-Darling normality test

data:  rstudent(modelo.dca)
A = 0.22365, p-value = 0.7973
lillie.test(rstudent(modelo.dca))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dca)
D = 0.11436, p-value = 0.7055
ks.test(rstudent(modelo.dca), "pnorm",
        alternative = "two.sided")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dca)
D = 0.12938, p-value = 0.8913
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dca))

    Cramer-von Mises normality test

data:  rstudent(modelo.dca)
W = 0.036876, p-value = 0.7224
pearson.test(rstudent(modelo.dca))

    Pearson chi-square normality test

data:  rstudent(modelo.dca)
P = 3.1, p-value = 0.5412
sf.test(rstudent(modelo.dca))

    Shapiro-Francia normality test

data:  rstudent(modelo.dca)
W = 0.98243, p-value = 0.9165

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento

ncvTest(modelo.dca)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.7649335, Df = 1, p = 0.38179
bptest(modelo.dca)

    studentized Breusch-Pagan test

data:  modelo.dca
BP = 1.887, df = 2, p-value = 0.3893
bptest(modelo.dca, studentize = F)

    Breusch-Pagan test

data:  modelo.dca
BP = 0.93752, df = 2, p-value = 0.6258
olsrr::ols_test_breusch_pagan(modelo.dca)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

                 Data                   
 ---------------------------------------
 Response : Rendimiento 
 Variables: fitted values of Rendimiento 

        Test Summary         
 ----------------------------
 DF            =    1 
 Chi2          =    0.7649335 
 Prob > Chi2   =    0.381789

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que no se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dca %>% gvlma()

Call:
lm(formula = Rendimiento ~ Variedad + Densidad, data = data)

Coefficients:
(Intercept)   Variedada2   Densidadb2  
       8.25         2.50        -3.50  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                      Value p-value                Decision
Global Stat        1.244611  0.8707 Assumptions acceptable.
Skewness           0.024166  0.8765 Assumptions acceptable.
Kurtosis           0.843892  0.3583 Assumptions acceptable.
Link Function      0.371747  0.5421 Assumptions acceptable.
Heteroscedasticity 0.004805  0.9447 Assumptions acceptable.

Análisis de varianza

\[Y_{ij} = \mu + \tau_{i} + \beta_{j} + \epsilon_{ij}\]

\[\hat{Y}_{ij} = \mu + \tau_{i} + \beta_{j}\]

Para el factor A (Variedad):

\(H_0: \tau_{A1} = \tau_{A2} = 0\)

\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)

Para el factor B (Densidad):

\(H_0: \beta_{B1} = \beta_{B2} = 0\)

\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)

anova(modelo.dca, test = "F")
Analysis of Variance Table

Response: Rendimiento
          Df Sum Sq Mean Sq F value   Pr(>F)   
Variedad   1  31.25  31.250  7.8996 0.012033 * 
Densidad   1  61.25  61.250 15.4833 0.001068 **
Residuals 17  67.25   3.956                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 1, 17)
[1] 4.451322

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 1, 17)
[1] 4.451322
shadeDist(qf(0.95, 1, 17), "df",1, 17, lower.tail = F)

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor A tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor A.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor B.

agricolae::cv.model(modelo.dca)
[1] 25.66374

Comparaciones de medias

Para los niveles del factor A:

  • A1 vs A2:

\(H_0: \mu_{A1} - \mu_{A2} = 0\)

\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)

Para los niveles del factor B:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

Prueba de LSD

Para los niveles del Factor A:

agricolae::LSD.test(modelo.dca, trt = "Variedad", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dca ~ "Variedad"

LSD t Test for Rendimiento 

Mean Square Error:  3.955882 

Variedad,  means and individual ( 95 %) CI

   Rendimiento      std  r      LCL       UCL Min Max
a1         6.5 2.273030 10 5.173014  7.826986   3  10
a2         9.0 3.018462 10 7.673014 10.326986   5  14

Alpha: 0.05 ; DF Error: 17
Critical Value of t: 2.109816 

least Significant Difference: 1.876641 

Treatments with the same letter are not significantly different.

   Rendimiento groups
a2         9.0      a
a1         6.5      b
agricolae::LSD.test(modelo.dca, trt = "Variedad", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dca ~ "Variedad"

LSD t Test for Rendimiento 

Mean Square Error:  3.955882 

Variedad,  means and individual ( 95 %) CI

   Rendimiento      std  r      LCL       UCL Min Max
a1         6.5 2.273030 10 5.173014  7.826986   3  10
a2         9.0 3.018462 10 7.673014 10.326986   5  14

Alpha: 0.05 ; DF Error: 17
Critical Value of t: 2.109816 

Comparison between treatments means

        difference pvalue signif.       LCL        UCL
a1 - a2       -2.5  0.012       * -4.376641 -0.6233591

Para los niveles del Factor B:

agricolae::LSD.test(modelo.dca, trt = "Densidad", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dca ~ "Densidad"

LSD t Test for Rendimiento 

Mean Square Error:  3.955882 

Densidad,  means and individual ( 95 %) CI

   Rendimiento      std  r      LCL       UCL Min Max
b1         9.5 2.592725 10 8.173014 10.826986   5  14
b2         6.0 2.054805 10 4.673014  7.326986   3  10

Alpha: 0.05 ; DF Error: 17
Critical Value of t: 2.109816 

least Significant Difference: 1.876641 

Treatments with the same letter are not significantly different.

   Rendimiento groups
b1         9.5      a
b2         6.0      b
agricolae::LSD.test(modelo.dca, trt = "Densidad", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dca ~ "Densidad"

LSD t Test for Rendimiento 

Mean Square Error:  3.955882 

Densidad,  means and individual ( 95 %) CI

   Rendimiento      std  r      LCL       UCL Min Max
b1         9.5 2.592725 10 8.173014 10.826986   5  14
b2         6.0 2.054805 10 4.673014  7.326986   3  10

Alpha: 0.05 ; DF Error: 17
Critical Value of t: 2.109816 

Comparison between treatments means

        difference pvalue signif.      LCL      UCL
b1 - b2        3.5 0.0011      ** 1.623359 5.376641

DCA factorial pxq con el paquete ExpDes

attach(data)
fat2.crd(factor1 = Variedad,
         factor2 =  Densidad,
         resp =  Rendimiento,
         quali = c(TRUE,
                   TRUE),
         mcomp = "tukey", 
         fac.names = c("Variedades",
                       "Densidad"),
         sigT = 0.05,
         sigF = 0.05,
         unfold = NULL)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Variedades 
FACTOR 2:  Densidad 
------------------------------------------------------------------------


Analysis of Variance Table
------------------------------------------------------------------------
                    DF     SS MS      Fc   Pr>Fc
Variedades           1  31.25  3  7.5758 0.01416
Densidad             1  61.25  5 14.8485 0.00141
Variedades*Densidad  1   1.25  2  0.3030 0.58959
Residuals           16  66.00  4                
Total               19 159.75  1                
------------------------------------------------------------------------
CV = 26.21 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.1811829 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------

No significant interaction: analyzing the simple effect
------------------------------------------------------------------------
Variedades
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    a2      9 
 b   a1      6.5 
------------------------------------------------------------------------

Densidad
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    b1      9.5 
 b   b2      6 
------------------------------------------------------------------------

DCA factorial pxq con tratamientos adicionales con el paquete ExpDes

data(ex8)
attach(ex8)
data(secaAd)
fat2.ad.crd(
  factor1 = inoculante,
  factor2 =  biodiesel,
  repet =  vaso,
  resp =  seca,
  respAd =  secaAd,
  quali = c(TRUE,FALSE),
  mcomp = "tukey",
  fac.names =
    c("Inoculant",
      "Biodiesel"), 
  sigT = 0.05, 
  sigF = 0.05)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Inoculant 
FACTOR 2:  Biodiesel 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                    DF      SS      MS      Fc  Pr>Fc
Inoculant            1 0.13054 0.13054 46.4363      0
Biodiesel            3 0.01575 0.00525  1.8671 0.1714
Inoculant*Biodiesel  3 0.02808 0.00936  3.3295 0.0429
Ad vs Factorial      1 0.07370  0.0737 26.2189  1e-04
Residuals           18 0.05060 0.00281               
Total               26 0.29867                       
------------------------------------------------------------------------
CV = 20.57 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.02718241 
WARNING: at 5% of significance, residuals can not be considered normal!
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
             Means  
Additional 0.11000 a
Factorial  0.27625 b
------------------------------------------------------------------------



Significant interaction: analyzing the interaction
------------------------------------------------------------------------

Analyzing  Inoculant  inside of each level of  Biodiesel 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                       DF      SS      MS      Fc  Pr>Fc
Biodiesel               3 0.01575 0.00525  1.8671 0.1714
Inoculant:Biodiesel 5   1 0.02042 0.02042  7.2628 0.0148
Inoculant:Biodiesel 10  1 0.09627 0.09627 34.2451      0
Inoculant:Biodiesel 15  1 0.03527 0.03527 12.5455 0.0023
Inoculant:Biodiesel 20  1 0.00667 0.00667  2.3715  0.141
Ad vs Factorial         1 0.07370  0.0737 26.2189  1e-04
Residuals              18 0.05060 0.00281               
Total                  26 0.29867                       
------------------------------------------------------------------------



 Inoculant  inside of the level  5  of  Biodiesel 
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    2   0.3133333 
 b   1   0.1966667 
------------------------------------------------------------------------


 Inoculant  inside of the level  10  of  Biodiesel 
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    2   0.4466667 
 b   1   0.1933333 
------------------------------------------------------------------------


 Inoculant  inside of the level  15  of  Biodiesel 
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    2   0.34 
 b   1   0.1866667 
------------------------------------------------------------------------


 Inoculant  inside of the level  20  of  Biodiesel 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels     Means
1      1 0.2333333
2      2 0.3000000
------------------------------------------------------------------------



Analyzing  Biodiesel  inside of each level of  Inoculant 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                            DF      SS      MS      Fc  Pr>Fc
Inoculant                    1 0.13054 0.13054 46.4363      0
Biodiesel:Inoculant esterco  3 0.00396 0.00132  0.4694 0.7073
Biodiesel:Inoculant mamona   3 0.03987 0.01329  4.7273 0.0133
Ad vs Factorial              1 0.07370  0.0737 26.2189  1e-04
Residuals                   18 0.05060 0.00281               
Total                       26 0.29867                       
------------------------------------------------------------------------



 Biodiesel  inside of the level  esterco  of  Inoculant 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels     Means
1     10 0.1933333
2     15 0.1866667
3     20 0.2333333
4      5 0.1966667
------------------------------------------------------------------------


 Biodiesel  inside of the level  mamona  of  Inoculant 
------------------------------------------------------------------------
Adjustment of polynomial models of regression
------------------------------------------------------------------------

Linear Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0  0.3867      0.0375     10.3137    0   
b1 -0.0029      0.0027     -1.0714 0.2982 
------------------------------------------

R2 of linear model
--------
0.080930
--------

Analysis of Variance of linear model
===========================================
              DF   SS     MS    Fc  p.value
-------------------------------------------
Linear Effect 1  0.0032 0.0032 1.15 0.29816
Lack of fit   2  0.0366 0.0183 6.52 0.00743
Residuals     18 0.0506 0.0028             
-------------------------------------------
------------------------------------------------------------------------

Quadratic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0  0.1700      0.0852     1.9949  0.0614 
b1  0.0404      0.0156     2.5983  0.0182 
b2 -0.0017      0.0006     -2.8312 0.0111 
------------------------------------------

R2 of quadratic model
--------
0.646100
--------

Analysis of Variance of quadratic model
==============================================
                 DF   SS     MS    Fc  p.value
----------------------------------------------
Linear Effect    1  0.0032 0.0032 1.15 0.29816
Quadratic Effect 1  0.0225 0.0225 8.02 0.01107
Lack of fit      1  0.0141 0.0141 5.02 0.03792
Residuals        18 0.0506 0.0028             
----------------------------------------------
------------------------------------------------------------------------

Cubic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 -0.3667      0.2543     -1.4420 0.1665 
b1  0.2111      0.0778     2.7144  0.0142 
b2 -0.0171      0.0069     -2.4834 0.0231 
b3  0.0004      0.0002     2.2401  0.0379 
------------------------------------------

R2 of cubic model
-
1
-

Analysis of Variance of cubic model
==============================================
                 DF   SS     MS    Fc  p.value
----------------------------------------------
Linear Effect    1  0.0032 0.0032 1.15 0.29816
Quadratic Effect 1  0.0225 0.0225 8.02 0.01107
Cubic Effect     1  0.0141 0.0141 5.02 0.03794
Lack of fit      0    0      0     0      1   
Residuals        18 0.0506 0.0028             
----------------------------------------------
------------------------------------------------------------------------
ggplot(ex8,
       aes(x = biodiesel,
           y = seca,
           colour = inoculante)) +
  geom_point() +
  geom_smooth()
ggplot(ex8 %>%
         group_by(biodiesel,
                  inoculante) %>%
         summarise(seca = mean(seca,
                               na.rm = T)),
       aes(x = biodiesel,
           y = seca,
           colour = inoculante)) +
  geom_point() +
  geom_smooth()
ggplot(ex8,
       aes(x = biodiesel,
           y = seca,
           colour = inoculante)) +
  geom_point() +
  geom_smooth(formula = y ~ x + I(x^2) + I(x^3),
              method = "lm")

Modelo final del efecto de Biodiesel dentro del Inoculante Mamona

\[ Y_i = \beta_0 + \beta_1*x_i + \beta_2*{x_i^2} + \beta_3*{x_i}^3+\epsilon_i\]

Dónde:

\(Y_i\) = Valor observado de seca.

\(x_i\) = Valor observado de Biodiesel dentro del Inoculante “Mamona”.

\(\beta_0\) = Intercepto = -0.3667.

\(\beta_1\) = Pendiente del termino lineal = -0.2111.

\(\beta_2\) = Pendiente del termino cuadrático = -0.0171.

\(\beta_3\) = Pendiente del termino cúbico = 0.0004.

Prueba de Dunnet

data <- ex8 %>% 
  dplyr::select(inoculante, biodiesel, vaso, seca) %>%
  dplyr::bind_rows(data.frame(inoculante = "ad",
                       biodiesel = 0,
                       vaso = 1:3,
                       seca = secaAd)) %>%
  dplyr::mutate(tratamiento =
                  interaction(inoculante,
                              biodiesel))
modelo <- lm(seca ~ tratamiento, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = seca ~ tratamiento, data = data)

Linear Hypotheses:
                       Estimate Std. Error t value Pr(>|t|)    
esterco.5 - ad.0 == 0   0.08667    0.04329   2.002  0.27562    
mamona.5 - ad.0 == 0    0.20333    0.04329   4.697  0.00121 ** 
esterco.10 - ad.0 == 0  0.08333    0.04329   1.925  0.31125    
mamona.10 - ad.0 == 0   0.33667    0.04329   7.777  < 0.001 ***
esterco.15 - ad.0 == 0  0.07667    0.04329   1.771  0.39150    
mamona.15 - ad.0 == 0   0.23000    0.04329   5.313  < 0.001 ***
esterco.20 - ad.0 == 0  0.12333    0.04329   2.849  0.05918 .  
mamona.20 - ad.0 == 0   0.19000    0.04329   4.389  0.00233 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Especificar el nivel del tratamiento adicional

data$tratamiento <-
  relevel(data$tratamiento,"ad.0")
modelo <- lm(seca ~ tratamiento, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = seca ~ tratamiento, data = data)

Linear Hypotheses:
                       Estimate Std. Error t value Pr(>|t|)    
esterco.5 - ad.0 == 0   0.08667    0.04329   2.002  0.27525    
mamona.5 - ad.0 == 0    0.20333    0.04329   4.697  0.00126 ** 
esterco.10 - ad.0 == 0  0.08333    0.04329   1.925  0.31122    
mamona.10 - ad.0 == 0   0.33667    0.04329   7.777  < 0.001 ***
esterco.15 - ad.0 == 0  0.07667    0.04329   1.771  0.39194    
mamona.15 - ad.0 == 0   0.23000    0.04329   5.313  < 0.001 ***
esterco.20 - ad.0 == 0  0.12333    0.04329   2.849  0.05939 .  
mamona.20 - ad.0 == 0   0.19000    0.04329   4.389  0.00236 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Diseño de bloques completos al azar en arreglo factorial con dos vías

Planeamiento

Crear un libro de campo con el paquete agricolae

trt <- c(3,4)

r <- 5

salida <- agricolae::design.ab(trt = trt,
                               r = r,
                               design = "rcbd",
                               serie = 3,
                               seed = 123,
                               kinds = "Super-Duper",
                               randomization = TRUE)
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Guardar el libro generado

write.table(salida$book,
            "books/rcbd2w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida$book,
           "books/rcbd2w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)

Guardar el libro de campo generado

write.table(salida %>% zigzag(),
            "books/rcbd2w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida %>% zigzag(),
           "books/rcbd2w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
# Función con errores, no debe replicar

agricolaeplotr::plot_design.factorial_rcbd(
  design = salida,
  factor_name = "A",
  reverse_y = TRUE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas")
# Función con errores, no debe replicar

agricolaeplotr::plot_design.factorial_rcbd(
  design = salida,
  factor_name = "B",
  reverse_y = TRUE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Columnas",
       y = "Filas")

Corrección de la función

plot_design.factorial_rcbd_corregido <- function (design, y = "row", factor_name = "A", width = 1, height = 1, 
  space_width = 0.95, space_height = 0.85, reverse_x = FALSE, 
  reverse_y = FALSE) 
{
  if (design$parameters$design == "factorial" && design$parameters$applied == 
    "rcbd") {
    test_string(y)
    test_input_reverse_x(reverse_x)
    test_input_reverse_y(reverse_y)
    test_input_height(height)
    test_input_width(width)
    test_input_height(space_height)
    test_input_width(space_width)
    test_name_in_column(factor_name, design)
    plots <- as.numeric(design$book[, 1])
    r <- length(levels(design$book$block))
    t <- length(design$book$plots)/r
    mat <- matrix(plots, byrow = TRUE, ncol = t)
    dims <- dim(mat)
    table <- as.table(mat)
    rownames(table) <- 1:dims[1]
    colnames(table) <- 1:dims[2]
    table <- as.data.frame(table)
    colnames(table) <- c("row", "col", "plots")
    table[, y] <- as.numeric(table[, y])
    table$col <- as.numeric(table$col)
    table <- merge(table, design$book, by.x = "plots", by.y = "plots")
    table$col <- rep(1:(length(table$plots)/max(as.numeric(table[, 
      y]))), times = max(as.numeric(table[, y])))
    table$col <- table$col * width
    table$row <- as.numeric(table[, y]) * height
    if (reverse_y == TRUE) {
      table$row <- abs(table$row - max(table$row)) + min(table$row)
    }
    if (reverse_x == TRUE) {
      table$col <- abs(table$col - max(table$col)) + min(table$col)
    }
    plt <- ggplot(table, aes_string(x = "col", y = y)) + 
      geom_tile(aes_string(fill = factor_name), width = width * 
        space_width, height = height * space_height) + 
      theme_bw() + theme(line = element_blank()) + geom_text(aes_string(label = "plots"), 
      colour = "black")
    return(plt)
  }
  else {
    stop(paste0("This is not the correct function for your experiment design.\n          A design from a factorial design with a\n                random complete block design (rcbd) is needed here.", 
      "you have a design of type ", design$parameters$design, 
      " and an applied type of ", design$parameters$applied, 
      "."))
  }
}
plot_design.factorial_rcbd_corregido(
  design = salida,
  factor_name = "A",
  reverse_y = TRUE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas")
# Función con errores, no debe replicar

plot_design.factorial_rcbd_corregido(
  design = salida,
  factor_name = "B",
  reverse_y = TRUE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Columnas",
       y = "Filas")

Crear un libro de campo con el paquete edibble

menu_factorial()
design("Factorial Design") %>%
  set_units(unit = 210) %>%
  set_trts(trt1 = 10,
           trt2 = 7) %>%
  allot_trts(~unit) %>%
  assign_trts("random", seed = 141) %>%
  serve_table()
rcbd <- takeout(menu_factorial(trt = c(3,4),
                              r = 5 ,
                              seed = 123))
rcbd %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
rcbd2 <- design("Factorial Design") %>%
  set_units(block = 5,
            unit = nested_in(block, 12)) %>%
  set_trts(trt1 = c("Amazon", "Coari x Lame", "Coari x Yangambí"),
           trt2 = c(0, 30, 60, 90)) %>%
  allot_trts(trt1 + trt2 ~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
rcbd2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
write.table(rcbd2 %>% as.data.frame(),
            "books/rcbd2w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(rcbd2 %>% as.data.frame(),
           "books/rcbd2w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
deggust::autoplot(rcbd2)
plot(rcbd2)

Análisis de DBCA en arreglo factorial con dos vías

Importación de datos

archivos <- list.files(pattern = "datos dbca 2w.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor) %>%
  mutate(Bloque = factor(Bloque))

Creación del modelo lineal

modelo.dbca <- lm(Respuesta ~ Variedad * Densidad + Bloque, data = data)

Definición del modelo

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*A_3 + \beta_3*B_2 + \beta_4*B_3 + \beta_5*Bloq_{II} + \beta_6*Bloq_{III} + \beta_7*Bloq_{IV} + \beta_8*Bloq_{V} + \beta_9*A_2*B_2 + \beta_{10}*A_3*B_2 + \beta_{11}*A_2*B_3 + \beta_{12}*A_3*B_3 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*A_3 + \beta_3*B_2 + \beta_4*B_3 + \beta_5*Bloq_{II} + \beta_6*Bloq_{III} + \beta_7*Bloq_{IV} + \beta_8*Bloq_{V} + \beta_9*A_2*B_2 + \beta_{10}*A_3*B_2 + \beta_{11}*A_2*B_3 + \beta_{12}*A_3*B_3\]

summary(modelo.dbca)

Call:
lm(formula = Respuesta ~ Variedad * Densidad + Bloque, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5111 -0.6000  0.2889  0.7111  2.4889 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)             8.7111     0.9171   9.499 7.85e-11 ***
Variedadv2              0.6000     1.0791   0.556  0.58206    
Variedadv3             -2.2000     1.0791  -2.039  0.04981 *  
Densidadd2              2.8000     1.0791   2.595  0.01417 *  
Densidadd3              0.6000     1.0791   0.556  0.58206    
Bloque2                 0.8889     0.8043   1.105  0.27733    
Bloque3                 2.0000     0.8043   2.487  0.01830 *  
Bloque4                 4.7778     0.8043   5.940 1.29e-06 ***
Bloque5                 3.7778     0.8043   4.697 4.79e-05 ***
Variedadv2:Densidadd2  -2.4000     1.5261  -1.573  0.12563    
Variedadv3:Densidadd2   1.6000     1.5261   1.048  0.30229    
Variedadv2:Densidadd3   0.6000     1.5261   0.393  0.69680    
Variedadv3:Densidadd3   5.4000     1.5261   3.539  0.00125 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.706 on 32 degrees of freedom
Multiple R-squared:  0.7388,    Adjusted R-squared:  0.6408 
F-statistic: 7.541 on 12 and 32 DF,  p-value: 2.346e-06

Verificación visual de los supuestos del modelo

performance::check_model(modelo.dbca)
ggResidpanel::resid_panel(modelo.dbca)
influence.measures(modelo.dbca)
Influence measures of
     lm(formula = Respuesta ~ Variedad * Densidad + Bloque, data = data) :

     dfb.1_  dfb.Vrd2  dfb.Vrd3  dfb.Dns2  dfb.Dns3  dfb.Blq2  dfb.Blq3
1   0.12604 -7.42e-02 -7.42e-02 -7.42e-02 -7.42e-02 -5.53e-02 -5.53e-02
2   0.10746 -1.03e-01 -1.03e-01 -1.03e-01 -1.03e-01  7.66e-02 -1.98e-17
3   0.07756 -7.42e-02 -7.42e-02 -7.42e-02 -7.42e-02  7.72e-17  5.53e-02
4  -0.70142  6.71e-01  6.71e-01  6.71e-01  6.71e-01 -8.00e-16 -4.08e-16
5   0.41263 -3.95e-01 -3.95e-01 -3.95e-01 -3.95e-01  2.13e-16 -6.31e-17
6  -0.20632 -2.24e-17 -2.63e-17 -3.95e-01 -8.12e-17  2.94e-01  2.94e-01
7  -0.02018  1.24e-17  2.51e-17  1.54e-01  6.03e-18  1.15e-01  7.74e-17
8  -0.05079  1.91e-16  3.19e-16  3.89e-01  1.94e-16 -1.35e-16  2.90e-01
9   0.00970  1.13e-17 -1.08e-17 -7.42e-02 -2.72e-19  2.38e-17  2.66e-18
10  0.00970 -2.33e-17 -4.17e-17 -7.42e-02 -1.56e-17  3.71e-17  2.64e-17
11  0.09275  1.65e-16  1.10e-16  1.40e-16  1.77e-01 -1.32e-01 -1.32e-01
12 -0.06190  1.35e-17  1.18e-17  4.16e-17  4.73e-01  3.53e-01  4.48e-17
13 -0.02319  1.07e-16  1.10e-16  3.11e-17  1.77e-01 -6.15e-18  1.32e-01
14  0.03685 -2.23e-16 -1.61e-16 -1.48e-16 -2.82e-01  7.57e-17 -1.01e-17
15  0.07249 -3.82e-16 -2.59e-16 -2.92e-16 -5.54e-01  1.20e-16  1.99e-17
16  0.09275  1.77e-01  9.12e-17  8.51e-17  1.11e-16 -1.32e-01 -1.32e-01
17  0.04069 -3.11e-01  2.23e-17  3.00e-17  2.37e-17 -2.32e-01 -7.60e-18
18 -0.02319  1.77e-01  8.04e-17  3.89e-17  9.34e-17 -8.61e-17  1.32e-01
19  0.07249 -5.54e-01 -2.74e-16 -2.91e-16 -2.51e-16  3.85e-16  7.95e-17
20 -0.06594  5.04e-01  1.89e-16  1.16e-16  1.09e-16 -2.10e-16  3.09e-17
21  0.31999  2.50e-16  2.66e-16  2.09e-16  2.15e-16 -4.56e-01 -4.56e-01
22  0.02018  9.64e-17  4.08e-17  4.31e-17  5.45e-17 -1.15e-01  1.16e-18
23  0.09644 -2.65e-17 -2.04e-16 -3.53e-17 -5.93e-18 -1.28e-16 -5.50e-01
24 -0.05158 -1.15e-17  9.24e-17  5.61e-17  1.59e-17 -1.02e-17  1.41e-17
25  0.01643 -5.66e-18 -2.72e-17 -1.48e-17  1.21e-17 -1.28e-17 -5.05e-18
26  0.35071  1.27e-16  1.93e-16  9.68e-17  2.68e-16 -5.00e-01 -5.00e-01
27  0.01343  6.58e-17  3.84e-17  3.83e-17  2.60e-17 -7.66e-02 -6.82e-19
28  0.13053  1.67e-16 -1.69e-16  2.23e-16 -1.46e-16 -6.93e-17 -7.44e-01
29 -0.05870 -9.97e-17  2.71e-17 -5.02e-17  5.56e-17 -1.10e-16 -6.44e-17
30  0.00970  2.11e-17 -4.25e-18  8.36e-18  1.96e-17 -3.09e-18 -1.10e-18
31 -0.06871 -1.15e-16 -1.31e-01 -3.93e-17 -6.68e-17  9.79e-02  9.79e-02
32  0.01343  1.80e-17 -1.03e-01  3.21e-17  2.93e-17 -7.66e-02  1.23e-17
33 -0.01643  3.31e-17  1.26e-01 -5.64e-17 -7.71e-18 -2.62e-17  9.36e-02
34 -0.02394  9.64e-17  1.83e-01  1.31e-17  6.66e-17 -5.88e-17 -6.56e-18
35  0.00970 -3.91e-17 -7.42e-02  2.13e-18 -9.42e-18  2.38e-17  1.17e-17
36 -0.12301 -1.35e-16 -1.17e-16 -2.09e-16 -9.05e-17  1.75e-01  1.75e-01
37 -0.00671 -2.09e-17 -6.07e-18  1.05e-17 -1.03e-17  3.83e-02  9.25e-18
38 -0.03685  1.49e-17  3.81e-17  4.51e-17  1.60e-17 -9.78e-18  2.10e-01
39 -0.01044  1.62e-17  1.35e-17  4.18e-17  1.14e-17 -1.17e-17  8.59e-18
40  0.02319 -2.67e-17 -3.20e-17 -6.41e-17  5.75e-18 -5.94e-18 -3.25e-17
41 -0.52213 -4.86e-16 -2.12e-16 -3.05e-16 -5.27e-16  7.44e-01  7.44e-01
42  0.01343  5.21e-17  4.25e-17  2.60e-17  1.84e-17 -7.66e-02  3.07e-18
43 -0.05079  7.92e-17  1.31e-16  2.03e-17  1.52e-16  4.04e-17  2.90e-01
44 -0.05870  9.05e-18  3.70e-17  5.86e-17  1.01e-16 -8.59e-17 -3.22e-17
45 -0.02394  2.06e-17  3.30e-17  9.63e-18  3.89e-17 -3.19e-17 -3.74e-17
    dfb.Blq4  dfb.Blq5 dfb.V2.D2 dfb.V3.D2 dfb.V2.D3 dfb.V3.D3   dffit cov.r
1  -5.53e-02 -5.53e-02  5.24e-02  5.24e-02  5.24e-02  5.24e-02  0.1260 2.090
2  -2.59e-17 -1.74e-17  7.26e-02  7.26e-02  7.26e-02  7.26e-02  0.1746 2.059
3   3.69e-17  4.93e-17  5.24e-02  5.24e-02  5.24e-02  5.24e-02  0.1260 2.090
4  -5.00e-01 -3.28e-16 -4.74e-01 -4.74e-01 -4.74e-01 -4.74e-01 -1.1398 0.593
5   1.50e-16  2.94e-01  2.79e-01  2.79e-01  2.79e-01  2.79e-01  0.6705 1.347
6   2.94e-01  2.94e-01  2.79e-01  2.79e-01  2.80e-17  2.60e-17 -0.6705 1.347
7   6.95e-17  8.03e-17 -1.09e-01 -1.09e-01 -1.24e-18  3.39e-18  0.2623 1.979
8  -1.47e-16 -8.78e-17 -2.75e-01 -2.75e-01 -1.66e-16 -2.13e-16  0.6603 1.365
9  -5.53e-02  1.03e-17  5.24e-02  5.24e-02 -5.26e-18 -4.88e-18 -0.1260 2.090
10  3.95e-17 -5.53e-02  5.24e-02  5.24e-02  1.92e-17  2.44e-17 -0.1260 2.090
11 -1.32e-01 -1.32e-01 -9.93e-17 -1.05e-16 -1.25e-01 -1.25e-01  0.3014 1.935
12 -1.72e-17  5.19e-18  6.59e-17  7.20e-18 -3.35e-01 -3.35e-01  0.8048 1.107
13 -6.96e-18 -5.84e-18 -5.39e-18 -1.08e-17 -1.25e-01 -1.25e-01  0.3014 1.935
14 -2.10e-01  3.40e-17  1.04e-16  6.00e-17  1.99e-01  1.99e-01 -0.4791 1.680
15  0.00e+00 -4.13e-01  2.26e-16  1.60e-16  3.92e-01  3.92e-01 -0.9424 0.876
16 -1.32e-01 -1.32e-01 -1.25e-01 -6.00e-17 -1.25e-01 -6.81e-17  0.3014 1.935
17  2.10e-17 -1.71e-17  2.20e-01 -2.06e-17  2.20e-01 -3.76e-17 -0.5290 1.597
18 -3.07e-17 -1.75e-17 -1.25e-01  6.50e-18 -1.25e-01 -5.64e-17  0.3014 1.935
19 -4.13e-01  1.03e-16  3.92e-01  1.54e-16  3.92e-01  1.89e-16 -0.9424 0.876
20  3.83e-17  3.76e-01 -3.57e-01 -1.00e-16 -3.57e-01 -7.19e-17  0.8572 1.017
21 -4.56e-01 -4.56e-01  4.33e-01 -3.18e-16 -1.79e-16 -1.95e-16  1.0400 0.728
22 -1.96e-17  2.37e-18 -1.09e-01 -9.21e-18 -6.64e-17 -3.89e-17 -0.2623 1.979
23 -5.34e-17  9.23e-17 -5.22e-01  8.44e-17  1.43e-17  5.66e-17 -1.2537 0.461
24  2.94e-01 -4.07e-17  2.79e-01 -9.13e-17  6.18e-18 -4.76e-17  0.6705 1.347
25 -2.87e-17 -9.36e-02 -8.88e-02  2.97e-17  5.83e-18  6.89e-18 -0.2135 2.027
26 -5.00e-01 -5.00e-01 -1.22e-16 -1.56e-16  4.74e-01 -2.28e-16  1.1398 0.593
27 -5.02e-18  5.18e-18 -3.51e-17 -2.96e-17 -7.26e-02 -2.14e-17 -0.1746 2.059
28  7.38e-17  1.82e-16 -2.20e-16 -1.76e-16 -7.06e-01  1.86e-16 -1.6969 0.146
29  3.35e-01 -1.01e-16  8.45e-17  5.07e-17  3.17e-01 -6.40e-17  0.7631 1.181
30 -1.13e-17 -5.53e-02 -2.87e-18  4.03e-18 -5.24e-02 -4.07e-18 -0.1260 2.090
31  9.79e-02  9.79e-02  5.19e-17  9.29e-02  5.99e-17  9.29e-02 -0.2233 2.018
32  1.05e-17  6.76e-19 -3.38e-17  7.26e-02 -3.90e-17  7.26e-02 -0.1746 2.059
33 -2.07e-17 -1.19e-17  4.14e-17 -8.88e-02  0.00e+00 -8.88e-02  0.2135 2.027
34  1.36e-01 -5.22e-18 -1.21e-17 -1.29e-01 -5.57e-17 -1.29e-01  0.3112 1.923
35  1.13e-17 -5.53e-02  1.46e-17  5.24e-02  2.25e-17  5.24e-02 -0.1260 2.090
36  1.75e-01  1.75e-01  1.55e-16 -1.66e-01  5.97e-17  5.68e-17 -0.3998 1.803
37  3.60e-18  5.63e-19  3.38e-18  3.63e-02  2.23e-17  5.63e-18  0.0872 2.108
38 -2.64e-17 -4.64e-17  0.00e+00  1.99e-01  2.26e-18 -6.18e-18  0.4791 1.680
39  5.95e-02 -6.13e-18 -2.63e-17  5.65e-02 -5.66e-18 -5.26e-18  0.1357 2.085
40 -4.04e-17 -1.32e-01  4.67e-17 -1.25e-01  5.32e-18  3.89e-18 -0.3014 1.935
41  7.44e-01  7.44e-01  3.31e-16  1.52e-16  3.79e-16 -7.06e-01 -1.6969 0.146
42 -4.71e-19  6.09e-18 -3.51e-17 -3.12e-17 -3.90e-17 -7.26e-02 -0.1746 2.059
43  2.51e-17 -5.97e-18  4.87e-18 -2.29e-32 -1.18e-16  2.75e-01  0.6603 1.365
44  3.35e-01 -6.60e-17  2.29e-18  3.41e-17  0.00e+00  3.17e-01  0.7631 1.181
45 -1.39e-17  1.36e-01 -1.39e-17 -5.39e-33 -2.78e-17  1.29e-01  0.3112 1.923
     cook.d   hat inf
1  0.001260 0.289    
2  0.002415 0.289    
3  0.001260 0.289    
4  0.093511 0.289    
5  0.034470 0.289    
6  0.034470 0.289    
7  0.005434 0.289    
8  0.033464 0.289    
9  0.001260 0.289    
10 0.001260 0.289    
11 0.007164 0.289    
12 0.048910 0.289    
13 0.007164 0.289    
14 0.017899 0.289    
15 0.065870 0.289    
16 0.007164 0.289    
17 0.021738 0.289    
18 0.007164 0.289    
19 0.065870 0.289    
20 0.055135 0.289    
21 0.079087 0.289    
22 0.005434 0.289    
23 0.110955 0.289    
24 0.034470 0.289    
25 0.003608 0.289    
26 0.093511 0.289    
27 0.002415 0.289    
28 0.186099 0.289    
29 0.044199 0.289    
30 0.001260 0.289    
31 0.003944 0.289    
32 0.002415 0.289    
33 0.003608 0.289    
34 0.007634 0.289    
35 0.001260 0.289    
36 0.012531 0.289    
37 0.000604 0.289    
38 0.017899 0.289    
39 0.001461 0.289    
40 0.007164 0.289    
41 0.186099 0.289    
42 0.002415 0.289    
43 0.033464 0.289    
44 0.044199 0.289    
45 0.007634 0.289    
influenceIndexPlot(modelo.dbca)

Cumplimiento de supuestos del modelo lineal general

Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dbca,
                 reps = 5000,
                 max.lag = 5)
 lag Autocorrelation D-W Statistic p-value
   1    -0.122529686      2.238735  0.6964
   2    -0.235172816      2.430874  0.4444
   3    -0.199332061      2.334500  0.3496
   4     0.002433206      1.862754  0.9584
   5     0.200567218      1.309637  0.0248
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dbca, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dbca
DW = 2.2387, p-value = 0.7084
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dbca))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dbca)
W = 0.95113, p-value = 0.05608
ad.test(rstudent(modelo.dbca))

    Anderson-Darling normality test

data:  rstudent(modelo.dbca)
A = 0.66817, p-value = 0.07577
lillie.test(rstudent(modelo.dbca))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dbca)
D = 0.10743, p-value = 0.2155
ks.test(rstudent(modelo.dbca), "pnorm",
        alternative = "two.sided")

    Exact one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dbca)
D = 0.096378, p-value = 0.7614
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))

    Cramer-von Mises normality test

data:  rstudent(modelo.dbca)
W = 0.10313, p-value = 0.09826
pearson.test(rstudent(modelo.dbca))

    Pearson chi-square normality test

data:  rstudent(modelo.dbca)
P = 14.333, p-value = 0.04556
sf.test(rstudent(modelo.dbca))

    Shapiro-Francia normality test

data:  rstudent(modelo.dbca)
W = 0.95355, p-value = 0.06612

Conclusión. A un nivel de significancia de 0.1, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento no es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento

ncvTest(modelo.dbca)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.4249817, Df = 1, p = 0.51446
bptest(modelo.dbca)

    studentized Breusch-Pagan test

data:  modelo.dbca
BP = 13.923, df = 12, p-value = 0.3057
bptest(modelo.dbca, studentize = F)

    Breusch-Pagan test

data:  modelo.dbca
BP = 13.778, df = 12, p-value = 0.3151
olsrr::ols_test_breusch_pagan(modelo.dbca)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

                Data                  
 -------------------------------------
 Response : Respuesta 
 Variables: fitted values of Respuesta 

        Test Summary         
 ----------------------------
 DF            =    1 
 Chi2          =    0.4249817 
 Prob > Chi2   =    0.5144617

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que no se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dbca %>% gvlma()

Call:
lm(formula = Respuesta ~ Variedad * Densidad + Bloque, data = data)

Coefficients:
          (Intercept)             Variedadv2             Variedadv3  
               8.7111                 0.6000                -2.2000  
           Densidadd2             Densidadd3                Bloque2  
               2.8000                 0.6000                 0.8889  
              Bloque3                Bloque4                Bloque5  
               2.0000                 4.7778                 3.7778  
Variedadv2:Densidadd2  Variedadv3:Densidadd2  Variedadv2:Densidadd3  
              -2.4000                 1.6000                 0.6000  
Variedadv3:Densidadd3  
               5.4000  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                       Value p-value                Decision
Global Stat        5.0988046  0.2773 Assumptions acceptable.
Skewness           2.7328520  0.0983 Assumptions acceptable.
Kurtosis           0.0008139  0.9772 Assumptions acceptable.
Link Function      2.1991572  0.1381 Assumptions acceptable.
Heteroscedasticity 0.1659815  0.6837 Assumptions acceptable.

Análisis de varianza

\[Y_{ijk} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + (\tau\beta)_{ij} + \epsilon_{ijk}\]

\[\hat{Y}_{ijk} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + (\tau\beta)_{ij}\] Dónde:

\(Y_{ijk}\) = Valor observado de la variable respuesta.

\(\hat{Y}_{ijk}\) = Valor ajustado de la variable respuesta.

\(\mu\) = Promedio observado de la variable respuesta.

\(\tau_{i}\) = Efecto del i-ésimo nivel del factor A.

\(\beta_{j}\) = Efecto del j-ésimo nivel del factor B.

\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor Bloque.

\((\tau\beta)_{ij}\) = Efecto de la interacción de los niveles del factor A dentro de un nivel del factor B o visceversa.

\(\epsilon_{ijk}\) = Residuo observado del modelo.

Pruebas de hipótesis

Para el factor A (Variedad):

\(H_0: \tau_{A1} = \tau_{A2} = \tau_{A3} = 0\)

\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)

Para el factor B (Densidad):

\(H_0: \beta_{B1} = \beta_{B2} = \beta_{B3} = 0\)

\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_{j} \neq 0\text{; en al menos un nivel del factor B.}\)

Para la interacción entre factor A y factor B:

\(H_0: (\tau\beta)_{A1B1} = (\tau\beta)_{A1B2} = (\tau\beta)_{A1B3} = (\tau\beta)_{A2B1} = (\tau\beta)_{A2B2} = (\tau\beta)_{A2B3} = (\tau\beta)_{A3B1} = (\tau\beta)_{A3B2} = (\tau\beta)_{A3B3} = 0\)

\(H_1: \text{En al menos una interacción entre el factor A y el factor B el } (\tau\beta) \text{ es diferente a los demás.}\)

\(H_1: (\tau\beta)_{ij} \neq 0\text{; en al menos una interacción entre el factor A y el factor B.}\)

anova(modelo.dbca, test = "F")
Analysis of Variance Table

Response: Respuesta
                  Df  Sum Sq Mean Sq F value    Pr(>F)    
Variedad           2   0.178   0.089  0.0305 0.9699553    
Densidad           2  65.911  32.956 11.3206 0.0001915 ***
Bloque             4 141.244  35.311 12.1298 4.121e-06 ***
Variedad:Densidad  4  56.089  14.022  4.8168 0.0037229 ** 
Residuals         32  93.156   2.911                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 2, 32)
[1] 3.294537

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 2, 32)
[1] 3.294537

Valor de la tabla de F para la interacción AxB con una significancia de 0.05.

qf(0.95, 4, 32)
[1] 2.668437

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor A tienen un efecto sobre el rendimiento estadísticamente similar a 0.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor B.

Con respecto a la interacción entre el Factor A y Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos una interacción entre un nivel del factor A y un nivel del factor B existe un efecto de antagonismo o sinergismo sobre el rendimiento.

agricolae::cv.model(modelo.dbca)
[1] 14.01075

Comparaciones de medias

Para los niveles del factor A:

  • A1 vs A2:

\(H_0: \mu_{A1} - \mu_{A2} = 0\)

\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)

  • A1 vs A3:

\(H_0: \mu_{A1} - \mu_{A3} = 0\)

\(H_1: \mu_{A1} - \mu_{A3} \neq 0\)

  • A2 vs A3:

\(H_0: \mu_{A2} - \mu_{A3} = 0\)

\(H_1: \mu_{A2} - \mu_{A3} \neq 0\)

Para los niveles del factor B:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

  • B1 vs B3:

\(H_0: \mu_{B1} - \mu_{B3} = 0\)

\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)

  • B2 vs B3:

\(H_0: \mu_{B2} - \mu_{B3} = 0\)

\(H_1: \mu_{B2} - \mu_{B3} \neq 0\)

Prueba de LSD

Para los niveles del Factor A:

agricolae::LSD.test(modelo.dbca, trt = "Variedad", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dbca ~ "Variedad"

LSD t Test for Respuesta 

Mean Square Error:  2.911111 

Variedad,  means and individual ( 95 %) CI

   Respuesta      std  r      LCL      UCL Min Max
v1  12.13333 2.133631 15 11.23599 13.03068   9  16
v2  12.13333 2.559762 15 11.23599 13.03068   9  17
v3  12.26667 3.788454 15 11.36932 13.16401   6  19

Alpha: 0.05 ; DF Error: 32
Critical Value of t: 2.036933 

least Significant Difference: 1.269041 

Treatments with the same letter are not significantly different.

   Respuesta groups
v3  12.26667      a
v1  12.13333      a
v2  12.13333      a
agricolae::LSD.test(modelo.dbca, trt = "Variedad", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dbca ~ "Variedad"

LSD t Test for Respuesta 

Mean Square Error:  2.911111 

Variedad,  means and individual ( 95 %) CI

   Respuesta      std  r      LCL      UCL Min Max
v1  12.13333 2.133631 15 11.23599 13.03068   9  16
v2  12.13333 2.559762 15 11.23599 13.03068   9  17
v3  12.26667 3.788454 15 11.36932 13.16401   6  19

Alpha: 0.05 ; DF Error: 32
Critical Value of t: 2.036933 

Comparison between treatments means

        difference pvalue signif.       LCL      UCL
v1 - v2  0.0000000 1.0000         -1.269041 1.269041
v1 - v3 -0.1333333 0.8319         -1.402374 1.135707
v2 - v3 -0.1333333 0.8319         -1.402374 1.135707

Para los niveles del Factor B:

agricolae::LSD.test(modelo.dbca, trt = "Densidad", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dbca ~ "Densidad"

LSD t Test for Respuesta 

Mean Square Error:  2.911111 

Densidad,  means and individual ( 95 %) CI

   Respuesta      std  r       LCL      UCL Min Max
d1  10.46667 2.386470 15  9.569319 11.36401   6  15
d2  13.00000 2.420153 15 12.102653 13.89735   9  16
d3  13.06667 3.034720 15 12.169319 13.96401   9  19

Alpha: 0.05 ; DF Error: 32
Critical Value of t: 2.036933 

least Significant Difference: 1.269041 

Treatments with the same letter are not significantly different.

   Respuesta groups
d3  13.06667      a
d2  13.00000      a
d1  10.46667      b
agricolae::LSD.test(modelo.dbca, trt = "Densidad", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dbca ~ "Densidad"

LSD t Test for Respuesta 

Mean Square Error:  2.911111 

Densidad,  means and individual ( 95 %) CI

   Respuesta      std  r       LCL      UCL Min Max
d1  10.46667 2.386470 15  9.569319 11.36401   6  15
d2  13.00000 2.420153 15 12.102653 13.89735   9  16
d3  13.06667 3.034720 15 12.169319 13.96401   9  19

Alpha: 0.05 ; DF Error: 32
Critical Value of t: 2.036933 

Comparison between treatments means

         difference pvalue signif.       LCL       UCL
d1 - d2 -2.53333333 0.0003     *** -3.802374 -1.264293
d1 - d3 -2.60000000 0.0002     *** -3.869041 -1.330959
d2 - d3 -0.06666667 0.9155         -1.335707  1.202374

Comparaciones de medias para las interacciones

Para los niveles del factor A dentro del nivel B1:

  • A1 vs A2:

\(H_0: \mu_{A1} - \mu_{A2} = 0\)

\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)

  • A1 vs A3:

\(H_0: \mu_{A1} - \mu_{A3} = 0\)

\(H_1: \mu_{A1} - \mu_{A3} \neq 0\)

  • A2 vs A3:

\(H_0: \mu_{A2} - \mu_{A3} = 0\)

\(H_1: \mu_{A2} - \mu_{A3} \neq 0\)

Para los niveles del factor B dentro del nivel A1:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

  • B1 vs B3:

\(H_0: \mu_{B1} - \mu_{B3} = 0\)

\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)

  • B2 vs B3:

\(H_0: \mu_{B2} - \mu_{B3} = 0\)

\(H_1: \mu_{B2} - \mu_{B3} \neq 0\)

NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.

Análisis de interacción entre variables cualitativas con el paquete ExpDes

attach(data)
fat2.rbd(factor1 = Variedad,
         factor2 =  Densidad,
         block = Bloque,
         resp =  Respuesta,
         quali = c(TRUE,
                   TRUE),
         mcomp = "lsd", 
         fac.names = c("Variedades",
                       "Densidad"),
         sigT = 0.05,
         sigF = 0.05)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Variedades 
FACTOR 2:  Densidad 
------------------------------------------------------------------------


Analysis of Variance Table
------------------------------------------------------------------------
                    DF     SS MS      Fc   Pr>Fc
Block                4 141.24  6 12.1298 0.00000
Variedades           2   0.18  2  0.0305 0.96996
Densidad             2  65.91  5 11.3206 0.00019
Variedades*Densidad  4  56.09  3  4.8168 0.00372
Residuals           32  93.16  4                
Total               44 356.58  1                
------------------------------------------------------------------------
CV = 14.01 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.09015993 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------



Significant interaction: analyzing the interaction
------------------------------------------------------------------------

Analyzing  Variedades  inside of each level of  Densidad 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                       DF        SS       MS      Fc  Pr.Fc
Block                   4 141.24444 35.31111 12.1298      0
Densidad                2  65.91111 32.95556 11.3206  2e-04
Variedades:Densidad d1  2  21.73333 10.86667  3.7328 0.0349
Variedades:Densidad d2  2   8.40000      4.2  1.4427 0.2512
Variedades:Densidad d3  2  26.13333 13.06667  4.4885 0.0191
Residuals              32  93.15556  2.91111               
Total                  44 356.57778                        
------------------------------------------------------------------------



 Variedades  inside of the level  d1  of  Densidad 
------------------------------------------------------------------------
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    2   11.6 
a    1   11 
 b   3   8.8 
------------------------------------------------------------------------


 Variedades  inside of the level  d2  of  Densidad 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
    Levels     Means
1        1      13.8
2        2      12.0
3        3      13.2
------------------------------------------------------------------------


 Variedades  inside of the level  d3  of  Densidad 
------------------------------------------------------------------------
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    3   14.8 
ab   2   12.8 
 b   1   11.6 
------------------------------------------------------------------------



Analyzing  Densidad  inside of each level of  Variedades 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                       DF        SS       MS      Fc  Pr.Fc
Block                   4 141.24444 35.31111 12.1298      0
Variedades              2   0.17778  0.08889  0.0305   0.97
Densidad:Variedades v1  2  21.73333 10.86667  3.7328 0.0349
Densidad:Variedades v2  2   3.73333  1.86667  0.6412 0.5333
Densidad:Variedades v3  2  96.53333 48.26667 16.5802      0
Residuals              32  93.15556  2.91111               
Total                  44 356.57778                        
------------------------------------------------------------------------



 Densidad  inside of the level  v1  of  Variedades 
------------------------------------------------------------------------
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    2   13.8 
 b   3   11.6 
 b   1   11 
------------------------------------------------------------------------


 Densidad  inside of the level  v2  of  Variedades 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
    Levels     Means
1        1      11.6
2        2      12.0
3        3      12.8
------------------------------------------------------------------------


 Densidad  inside of the level  v3  of  Variedades 
------------------------------------------------------------------------
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    3   14.8 
a    2   13.2 
 b   1   8.8 
------------------------------------------------------------------------

DBCA factorial pxq con tratamientos adicionales con el paquete ExpDes

data(ex7)
attach(ex7)
data(est21Ad)
fat2.ad.rbd(factor1 = periodo,
            factor2 = nivel,
            block = bloco,
            resp = est21,
            respAd =  est21Ad,
            quali = c(TRUE, FALSE),
            mcomp = "lsd",
            fac.names = c("Period", "Level"), 
            sigT = 0.05, 
            sigF = 0.05,
            unfold = NULL)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Period 
FACTOR 2:  Level 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                DF        SS        MS      Fc  Pr>Fc
Block            3   18.4634   6.15447  0.8717 0.4609
Period           4  981.2159 245.30397 34.7426      0
Level            3  394.8364 131.61212 18.6403      0
Period*Level    12  174.7855  14.56546  2.0629 0.0336
Ad vs Factorial  1  270.3703 270.37029 38.2928      0
Residuals       60  423.6366   7.06061               
Total           83 2263.3081                         
------------------------------------------------------------------------
CV = 11.48 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.1854518 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
             Means  
Additional 31.1700 a
Factorial  22.7455 b
------------------------------------------------------------------------
fat2.ad.rbd(factor1 = periodo,
            factor2 = nivel,
            block = bloco,
            resp = est21,
            respAd =  est21Ad,
            quali = c(TRUE, FALSE),
            mcomp = "lsd",
            fac.names = c("Period", "Level"), 
            sigT = 0.05, 
            sigF = 0.05,
            unfold = 0)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Period 
FACTOR 2:  Level 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                DF        SS        MS      Fc  Pr>Fc
Block            3   18.4634   6.15447  0.8717 0.4609
Period           4  981.2159 245.30397 34.7426      0
Level            3  394.8364 131.61212 18.6403      0
Period*Level    12  174.7855  14.56546  2.0629 0.0336
Ad vs Factorial  1  270.3703 270.37029 38.2928      0
Residuals       60  423.6366   7.06061               
Total           83 2263.3081                         
------------------------------------------------------------------------
CV = 11.48 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.1854518 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
             Means  
Additional 31.1700 a
Factorial  22.7455 b
------------------------------------------------------------------------
fat2.ad.rbd(factor1 = periodo,
            factor2 = nivel,
            block = bloco,
            resp = est21,
            respAd =  est21Ad,
            quali = c(TRUE, FALSE),
            mcomp = "lsd",
            fac.names = c("Period", "Level"), 
            sigT = 0.05, 
            sigF = 0.05,
            unfold = 1)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Period 
FACTOR 2:  Level 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                DF        SS        MS      Fc  Pr>Fc
Block            3   18.4634   6.15447  0.8717 0.4609
Period           4  981.2159 245.30397 34.7426      0
Level            3  394.8364 131.61212 18.6403      0
Period*Level    12  174.7855  14.56546  2.0629 0.0336
Ad vs Factorial  1  270.3703 270.37029 38.2928      0
Residuals       60  423.6366   7.06061               
Total           83 2263.3081                         
------------------------------------------------------------------------
CV = 11.48 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.1854518 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
             Means  
Additional 31.1700 a
Factorial  22.7455 b
------------------------------------------------------------------------

No significant interaction: analyzing the simple effect
------------------------------------------------------------------------
Period
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    7-14DAE     26.52563 
a    7-21DAE     25.73813 
 b   0-7DAE      23.4275 
  c      0-14DAE     21.235 
   d     0-21DAE     16.80125 
------------------------------------------------------------------------

Level
Adjustment of polynomial models of regression
------------------------------------------------------------------------

Linear Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 27.6000      0.7277     37.9278    0   
b1 -0.9709      0.1329     -7.3078    0   
------------------------------------------

R2 of linear model
--------
0.954975
--------

Analysis of Variance of linear model
===============================================
              DF    SS       MS     Fc  p.value
-----------------------------------------------
Linear Effect 1  377.0587 377.0587 53.4    0   
Lack of fit   2  17.7777   8.8888  1.26 0.29135
Residuals     60 423.6366  7.0606              
-----------------------------------------------
------------------------------------------------------------------------

Quadratic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 25.2737      1.6541     15.2796    0   
b1  0.1922      0.7545     0.2548  0.7998 
b2 -0.1163      0.0743     -1.5661 0.1226 
------------------------------------------

R2 of quadratic model
--------
0.998832
--------

Analysis of Variance of quadratic model
==================================================
                 DF    SS       MS     Fc  p.value
--------------------------------------------------
Linear Effect    1  377.0587 377.0587 53.4    0   
Quadratic Effect 1  17.3166  17.3166  2.45 0.12259
Lack of fit      1   0.4610   0.4610  0.07 0.79918
Residuals        60 423.6366  7.0606              
--------------------------------------------------
------------------------------------------------------------------------

Cubic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 24.0855      4.9355     4.8801  0.00001
b1  1.1372      3.7741     0.3013  0.7642 
b2 -0.3285      0.8337     -0.3940 0.6950 
b3  0.0141      0.0554     0.2555  0.7992 
------------------------------------------

R2 of cubic model
-
1
-

Analysis of Variance of cubic model
==================================================
                 DF    SS       MS     Fc  p.value
--------------------------------------------------
Linear Effect    1  377.0587 377.0587 53.4    0   
Quadratic Effect 1  17.3166  17.3166  2.45 0.12259
Cubic Effect     1   0.4610   0.4610  0.07 0.79918
Lack of fit      0     0        0      0      1   
Residuals        60 423.6366  7.0606              
--------------------------------------------------
------------------------------------------------------------------------




Significant interaction: analyzing the interaction
------------------------------------------------------------------------

Analyzing  Period  inside of each level of  Level 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                DF         SS        MS      Fc  Pr>Fc
Block            3   18.46340   6.15447  0.8717 0.4609
Level            3  394.83637 131.61212 18.6403      0
Period:Level 2   4   68.46548  17.11637  2.4242 0.0578
Period:Level 4   4  219.04878  54.76219   7.756      0
Period:Level 6   4  246.36947  61.59237  8.7234      0
Period:Level 8   4  622.11768 155.52942 22.0278      0
Ad vs Factorial  1  270.37029 270.37029 38.2928      0
Residuals       60  423.63660   7.06061               
Total           83 2263.30807                         
------------------------------------------------------------------------



 Period  inside of the level  2  of  Level 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
    Levels     Means
1  0-14DAE    25.180
2  0-21DAE    22.165
3   0-7DAE    24.695
4  7-14DAE    27.885
5  7-21DAE    25.870
------------------------------------------------------------------------


 Period  inside of the level  4  of  Level 
------------------------------------------------------------------------
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    7-21DAE     27.93 
a    7-14DAE     27.015 
ab   0-7DAE      24.835 
 b   0-14DAE     23.02 
  c      0-21DAE     18.6175 
------------------------------------------------------------------------


 Period  inside of the level  6  of  Level 
------------------------------------------------------------------------
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    7-14DAE     26.21 
a    7-21DAE     25.2775 
ab   0-7DAE      22.63 
 bc      0-14DAE     19.905 
  c      0-21DAE     16.6675 
------------------------------------------------------------------------


 Period  inside of the level  8  of  Level 
------------------------------------------------------------------------
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    7-14DAE     24.9925 
a    7-21DAE     23.875 
a    0-7DAE      21.55 
 b   0-14DAE     16.835 
  c      0-21DAE     9.755 
------------------------------------------------------------------------



Analyzing  Level  inside of each level of  Period 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                     DF         SS        MS      Fc  Pr>Fc
Block                 3   18.46340   6.15447  0.8717 0.4609
Period                4  981.21587 245.30397 34.7426      0
Level:Period 0-14DAE  3  159.51260  53.17087  7.5306  2e-04
Level:Period 0-21DAE  3  326.94442 108.98147 15.4351      0
Level:Period 0-7DAE   3   30.99450   10.3315  1.4633 0.2336
Level:Period 7-14DAE  3   18.14992   6.04997  0.8569 0.4685
Level:Period 7-21DAE  3   34.02047  11.34016  1.6061 0.1974
Ad vs Factorial       1  270.37029 270.37029 38.2928      0
Residuals            60  423.63660   7.06061               
Total                83 2263.30807                         
------------------------------------------------------------------------



 Level  inside of the level  0-14DAE  of  Period 
------------------------------------------------------------------------
Adjustment of polynomial models of regression
------------------------------------------------------------------------

Linear Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 28.2725      1.6272     17.3751    0   
b1 -1.4075      0.2971     -4.7378 0.00001
------------------------------------------

R2 of linear model
--------
0.993555
--------

Analysis of Variance of linear model
================================================
              DF    SS       MS     Fc   p.value
------------------------------------------------
Linear Effect 1  158.4845 158.4845 22.45  1e-05 
Lack of fit   2   1.0281   0.5140  0.07  0.92986
Residuals     60 423.6366  7.0606               
------------------------------------------------
------------------------------------------------------------------------

Quadratic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 27.1350      3.6986     7.3365     0   
b1 -0.8388      1.6871     -0.4972 0.6209 
b2 -0.0569      0.1661     -0.3425 0.7332 
------------------------------------------

R2 of quadratic model
--------
0.998746
--------

Analysis of Variance of quadratic model
===================================================
                 DF    SS       MS     Fc   p.value
---------------------------------------------------
Linear Effect    1  158.4845 158.4845 22.45  1e-05 
Quadratic Effect 1   0.8281   0.8281  0.12  0.7332 
Lack of fit      1   0.2000   0.2000  0.03  0.86691
Residuals        60 423.6366  7.0606               
---------------------------------------------------
------------------------------------------------------------------------

Cubic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 25.3850     11.0361     2.3002  0.0249 
b1  0.5529      8.4391     0.0655  0.9480 
b2 -0.3694      1.8642     -0.1981 0.8436 
b3  0.0208      0.1238     0.1683  0.8669 
------------------------------------------

R2 of cubic model
-
1
-

Analysis of Variance of cubic model
===================================================
                 DF    SS       MS     Fc   p.value
---------------------------------------------------
Linear Effect    1  158.4845 158.4845 22.45  1e-05 
Quadratic Effect 1   0.8281   0.8281  0.12  0.7332 
Cubic Effect     1   0.2000   0.2000  0.03  0.86691
Lack of fit      0     0        0       0      1   
Residuals        60 423.6366  7.0606               
---------------------------------------------------
------------------------------------------------------------------------


 Level  inside of the level  0-21DAE  of  Period 
------------------------------------------------------------------------
Adjustment of polynomial models of regression
------------------------------------------------------------------------

Linear Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 26.5963      1.6272     16.3450    0   
b1 -1.9590      0.2971     -6.5941    0   
------------------------------------------

R2 of linear model
--------
0.939042
--------

Analysis of Variance of linear model
================================================
              DF    SS       MS     Fc   p.value
------------------------------------------------
Linear Effect 1  307.0145 307.0145 43.48    0   
Lack of fit   2  19.9299   9.9650  1.41  0.25179
Residuals     60 423.6366  7.0606               
------------------------------------------------
------------------------------------------------------------------------

Quadratic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 22.3900      3.6986     6.0536     0   
b1  0.1441      1.6871     0.0854  0.9322 
b2 -0.2103      0.1661     -1.2664 0.2103 
------------------------------------------

R2 of quadratic model
--------
0.973675
--------

Analysis of Variance of quadratic model
===================================================
                 DF    SS       MS     Fc   p.value
---------------------------------------------------
Linear Effect    1  307.0145 307.0145 43.48    0   
Quadratic Effect 1  11.3232  11.3232   1.6  0.21027
Lack of fit      1   8.6067   8.6067  1.22  0.27397
Residuals        60 423.6366  7.0606               
---------------------------------------------------
------------------------------------------------------------------------

Cubic Model
==========================================
   Estimate Standard.Error   tc    p.value
------------------------------------------
b0 33.8700     11.0361     3.0690  0.0032 
b1 -8.9852      8.4391     -1.0647 0.2913 
b2  1.8397      1.8642     0.9869  0.3277 
b3 -0.1367      0.1238     -1.1041 0.2740 
------------------------------------------

R2 of cubic model
-
1
-

Analysis of Variance of cubic model
===================================================
                 DF    SS       MS     Fc   p.value
---------------------------------------------------
Linear Effect    1  307.0145 307.0145 43.48    0   
Quadratic Effect 1  11.3232  11.3232   1.6  0.21027
Cubic Effect     1   8.6067   8.6067  1.22  0.27397
Lack of fit      0     0        0       0      1   
Residuals        60 423.6366  7.0606               
---------------------------------------------------
------------------------------------------------------------------------


 Level  inside of the level  0-7DAE  of  Period 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
    Levels     Means
1        2    24.695
2        4    24.835
3        6    22.630
4        8    21.550
------------------------------------------------------------------------


 Level  inside of the level  7-14DAE  of  Period 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
    Levels     Means
1        2   27.8850
2        4   27.0150
3        6   26.2100
4        8   24.9925
------------------------------------------------------------------------


 Level  inside of the level  7-21DAE  of  Period 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
    Levels     Means
1        2   25.8700
2        4   27.9300
3        6   25.2775
4        8   23.8750
------------------------------------------------------------------------
fat2.ad.rbd(factor1 = periodo,
            factor2 = nivel,
            block = bloco,
            resp = est21,
            respAd =  est21Ad,
            quali = c(TRUE, FALSE),
            mcomp = "lsd",
            fac.names = c("Period", "Level"), 
            sigT = 0.05, 
            sigF = 0.05,
            unfold = 2)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Period 
FACTOR 2:  Level 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                DF        SS        MS      Fc  Pr>Fc
Block            3   18.4634   6.15447  0.8717 0.4609
Period           4  981.2159 245.30397 34.7426      0
Level            3  394.8364 131.61212 18.6403      0
Period*Level    12  174.7855  14.56546  2.0629 0.0336
Ad vs Factorial  1  270.3703 270.37029 38.2928      0
Residuals       60  423.6366   7.06061               
Total           83 2263.3081                         
------------------------------------------------------------------------
CV = 11.48 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.1854518 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
             Means  
Additional 31.1700 a
Factorial  22.7455 b
------------------------------------------------------------------------

Prueba de Dunnet

data <- ex7 %>% 
  dplyr::select(periodo, nivel, bloco, est21) %>%
  dplyr::bind_rows(data.frame(periodo = "0",
                       nivel = 0,
                       bloco = 1:4,
                       est21 = est21Ad)) %>%
  dplyr::mutate(tratamiento =
                  interaction(periodo,
                              nivel))
modelo <- lm(est21 ~ tratamiento + bloco, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = est21 ~ tratamiento + bloco, data = data)

Linear Hypotheses:
                     Estimate Std. Error t value Pr(>|t|)    
0-14DAE.2 - 0.0 == 0   -5.990      1.851  -3.236   0.0273 *  
0-21DAE.2 - 0.0 == 0   -9.005      1.851  -4.864    <0.01 ***
0-7DAE.2 - 0.0 == 0    -6.475      1.851  -3.498   0.0129 *  
7-14DAE.2 - 0.0 == 0   -3.285      1.851  -1.774   0.5671    
7-21DAE.2 - 0.0 == 0   -5.300      1.851  -2.863   0.0704 .  
0-14DAE.4 - 0.0 == 0   -8.150      1.851  -4.402    <0.01 ***
0-21DAE.4 - 0.0 == 0  -12.553      1.851  -6.780    <0.01 ***
0-7DAE.4 - 0.0 == 0    -6.335      1.851  -3.422   0.0158 *  
7-14DAE.4 - 0.0 == 0   -4.155      1.851  -2.244   0.2659    
7-21DAE.4 - 0.0 == 0   -3.240      1.851  -1.750   0.5859    
0-14DAE.6 - 0.0 == 0  -11.265      1.851  -6.085    <0.01 ***
0-21DAE.6 - 0.0 == 0  -14.503      1.851  -7.834    <0.01 ***
0-7DAE.6 - 0.0 == 0    -8.540      1.851  -4.613    <0.01 ***
7-14DAE.6 - 0.0 == 0   -4.960      1.851  -2.679   0.1082    
7-21DAE.6 - 0.0 == 0   -5.892      1.851  -3.183   0.0314 *  
0-14DAE.8 - 0.0 == 0  -14.335      1.851  -7.743    <0.01 ***
0-21DAE.8 - 0.0 == 0  -21.415      1.851 -11.567    <0.01 ***
0-7DAE.8 - 0.0 == 0    -9.620      1.851  -5.196    <0.01 ***
7-14DAE.8 - 0.0 == 0   -6.177      1.851  -3.337   0.0204 *  
7-21DAE.8 - 0.0 == 0   -7.295      1.851  -3.940    <0.01 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Especificar el nivel del tratamiento adicional

data$tratamiento <-
  relevel(data$tratamiento,"0.0")
modelo <- lm(est21 ~ tratamiento + bloco, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = est21 ~ tratamiento + bloco, data = data)

Linear Hypotheses:
                     Estimate Std. Error t value Pr(>|t|)    
0-14DAE.2 - 0.0 == 0   -5.990      1.851  -3.236  0.02714 *  
0-21DAE.2 - 0.0 == 0   -9.005      1.851  -4.864  < 0.001 ***
0-7DAE.2 - 0.0 == 0    -6.475      1.851  -3.498  0.01299 *  
7-14DAE.2 - 0.0 == 0   -3.285      1.851  -1.774  0.56748    
7-21DAE.2 - 0.0 == 0   -5.300      1.851  -2.863  0.07023 .  
0-14DAE.4 - 0.0 == 0   -8.150      1.851  -4.402  < 0.001 ***
0-21DAE.4 - 0.0 == 0  -12.553      1.851  -6.780  < 0.001 ***
0-7DAE.4 - 0.0 == 0    -6.335      1.851  -3.422  0.01625 *  
7-14DAE.4 - 0.0 == 0   -4.155      1.851  -2.244  0.26613    
7-21DAE.4 - 0.0 == 0   -3.240      1.851  -1.750  0.58575    
0-14DAE.6 - 0.0 == 0  -11.265      1.851  -6.085  < 0.001 ***
0-21DAE.6 - 0.0 == 0  -14.503      1.851  -7.834  < 0.001 ***
0-7DAE.6 - 0.0 == 0    -8.540      1.851  -4.613  < 0.001 ***
7-14DAE.6 - 0.0 == 0   -4.960      1.851  -2.679  0.10827    
7-21DAE.6 - 0.0 == 0   -5.892      1.851  -3.183  0.03118 *  
0-14DAE.8 - 0.0 == 0  -14.335      1.851  -7.743  < 0.001 ***
0-21DAE.8 - 0.0 == 0  -21.415      1.851 -11.567  < 0.001 ***
0-7DAE.8 - 0.0 == 0    -9.620      1.851  -5.196  < 0.001 ***
7-14DAE.8 - 0.0 == 0   -6.177      1.851  -3.337  0.02046 *  
7-21DAE.8 - 0.0 == 0   -7.295      1.851  -3.940  0.00337 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Diseño cuadrado latino en arreglo factorial con dos vías

Planeamiento

Crear un libro de campo con el paquete agricolae

trt <- c(3,4)

salida <- agricolae::design.ab(trt = trt,
                               design = "lsd",
                               serie = 3,
                               seed = 123,
                               kinds = "Super-Duper",
                               randomization = TRUE)
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Guardar el libro generado

write.table(salida$book,
            "books/lsd2w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida$book,
           "books/lsd2w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)

Guardar el libro de campo generado

write.table(salida %>% zigzag(),
            "books/lsd2w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida %>% zigzag(),
           "books/lsd2w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
agricolaeplotr::plot_design.factorial_lsd(
  design = salida,
  factor_name = "A",
  reverse_y = TRUE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas")
agricolaeplotr::plot_design.factorial_lsd(
  design = salida,
  factor_name = "B",
  reverse_y = TRUE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Columnas",
       y = "Filas")

Crear un libro de campo con el paquete edibble

menu_factorial()
design("Factorial Design") %>%
  set_units(unit = 540) %>%
  set_trts(trt1 = 10,
           trt2 = 9) %>%
  allot_trts(~unit) %>%
  assign_trts("random", seed = 579) %>%
  serve_table()
lsd <- takeout(menu_factorial(trt = c(3,4),
                              r = 12,
                              seed = 123))
lsd %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
lsd2 <- design("Factorial Design") %>%
  set_units(row = 12,
            col = 12,
            unit = crossed_by(row, col)) %>%
  set_trts(trt1 = c("Amazon", "Coari x Lame", "Coari x Yangambí"),
           trt2 = c(0, 30, 60, 90)) %>%
  allot_trts(trt1 + trt2 ~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
lsd2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
write.table(lsd2 %>% as.data.frame(),
            "books/rcbd2w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(lsd2 %>% as.data.frame(),
           "books/rcbd2w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
deggust::autoplot(lsd2)
plot(lsd2)

Análisis de DCL en arreglo factorial con dos vías

Importación de datos

archivos <- list.files(pattern = "datos dcl 2w.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor) %>%
  mutate(row = factor(row),
         col = factor(col))

Creación del modelo lineal con interacción

modelo.dcl <- lm(rdto ~ sp * trt + row + col, data = data)

Definición del modelo

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*F_{2} + \beta_4*F_{3} + \beta_5*F_{4} + \beta_6*C_2 + \beta_{7}*C_3 + \beta_{8}*C_4 + \beta_{9}*A_2*B_2 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*F_{2} + \beta_4*F_{3} + \beta_5*F_{4} + \beta_6*C_2 + \beta_{7}*C_3 + \beta_{8}*C_4 + \beta_{9}*A_2*B_2\]

summary(modelo.dcl)

Call:
lm(formula = rdto ~ sp * trt + row + col, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-1.250 -0.425  0.050  0.375  1.375 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.2000     1.0840   7.565 0.000277 ***
spsp2         0.6000     0.9695   0.619 0.558761    
trtB         -3.9250     0.9695  -4.048 0.006739 ** 
row2          2.1000     0.9695   2.166 0.073468 .  
row3          1.3500     0.9695   1.392 0.213209    
row4          0.8000     0.9695   0.825 0.440858    
col2          0.6500     0.9695   0.670 0.527535    
col3         -0.3500     0.9695  -0.361 0.730471    
col4          0.7500     0.9695   0.774 0.468560    
spsp2:trtB    0.6000     1.3711   0.438 0.676992    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.371 on 6 degrees of freedom
Multiple R-squared:  0.8593,    Adjusted R-squared:  0.6484 
F-statistic: 4.073 on 9 and 6 DF,  p-value: 0.05071

Creación del modelo lineal sin interacción

modelo.dcl <- lm(rdto ~ sp + trt + row + col, data = data)

Definición del modelo

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*F_{2} + \beta_4*F_{3} + \beta_5*F_{4} + \beta_6*C_2 + \beta_{7}*C_3 + \beta_{8}*C_4 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*F_{2} + \beta_4*F_{3} + \beta_5*F_{4} + \beta_6*C_2 + \beta_{7}*C_3 + \beta_{8}*C_4\]

summary(modelo.dcl)

Call:
lm(formula = rdto ~ sp + trt + row + col, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4000 -0.4813 -0.0500  0.4813  1.4000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.0500     0.9671   8.324 7.07e-05 ***
spsp2         0.9000     0.6448   1.396 0.205429    
trtB         -3.6250     0.6448  -5.622 0.000797 ***
row2          2.1000     0.9118   2.303 0.054743 .  
row3          1.3500     0.9118   1.481 0.182274    
row4          0.8000     0.9118   0.877 0.409391    
col2          0.6500     0.9118   0.713 0.499003    
col3         -0.3500     0.9118  -0.384 0.712486    
col4          0.7500     0.9118   0.823 0.437896    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.29 on 7 degrees of freedom
Multiple R-squared:  0.8549,    Adjusted R-squared:  0.689 
F-statistic: 5.154 on 8 and 7 DF,  p-value: 0.02186

Verificación visual de los supuestos del modelo

performance::check_model(modelo.dcl)
ggResidpanel::resid_panel(modelo.dcl)
influence.measures(modelo.dcl)
Influence measures of
     lm(formula = rdto ~ sp + trt + row + col, data = data) :

    dfb.1_ dfb.sps2 dfb.trtB  dfb.row2  dfb.row3  dfb.row4  dfb.col2  dfb.col3
1  -0.0616   0.0205   0.0205  2.90e-02  2.90e-02  2.90e-02  2.90e-02  2.90e-02
2   0.5984  -0.5984   0.5984 -8.46e-01 -8.46e-01 -8.46e-01  8.46e-01 -1.20e-16
3  -0.0411  -0.0411   0.0411  5.81e-02  5.81e-02  5.81e-02  3.35e-17 -5.81e-02
4  -0.1672  -0.5016  -0.5016  7.09e-01  7.09e-01  7.09e-01 -1.74e-16 -3.86e-16
5  -0.0205  -0.0205   0.0205 -2.90e-02  2.42e-18  5.55e-18  2.90e-02  2.90e-02
6  -0.5324   0.5324   0.5324  7.53e-01 -7.33e-17 -4.73e-33  7.53e-01 -1.81e-16
7   0.0275  -0.0824  -0.0824  1.16e-01  1.79e-17  1.20e-17  1.49e-17  1.16e-01
8   0.2113   0.6340  -0.6340 -8.97e-01  2.61e-17 -8.80e-18  2.40e-16  5.49e-16
9   0.0449   0.1347   0.1347  2.44e-17  1.91e-01 -6.54e-18 -1.91e-01 -1.91e-01
10 -0.1432   0.4297   0.4297  1.56e-16 -6.08e-01  9.54e-17 -6.08e-01 -3.20e-18
11  0.0449   0.1347  -0.1347  4.89e-17 -1.91e-01  3.27e-17  1.22e-17 -1.91e-01
12 -0.1432   0.4297  -0.4297  3.90e-17  6.08e-01  1.25e-16 -2.64e-16 -2.07e-16
13 -0.0928   0.0928  -0.0928  5.95e-18  1.78e-17 -1.31e-01  1.31e-01  1.31e-01
14  0.2441  -0.7323   0.7323  1.37e-16  2.11e-16 -1.04e+00 -1.04e+00 -2.98e-17
15 -0.0928   0.0928   0.0928 -1.44e-17 -1.78e-17  1.31e-01  8.41e-17  1.31e-01
16  0.2441  -0.7323  -0.7323 -8.99e-17 -7.04e-17  1.04e+00 -4.69e-17 -1.41e-16
    dfb.col4   dffit cov.r   cook.d   hat inf
1   2.90e-02 -0.0616 9.112 0.000491 0.563   *
2  -6.64e-17  1.7952 0.395 0.294674 0.562    
3   2.51e-17 -0.1232 8.992 0.001964 0.562   *
4  -7.09e-01 -1.5048 0.903 0.226927 0.562    
5   2.90e-02 -0.0616 9.112 0.000491 0.562   *
6  -2.96e-17  1.5973 0.699 0.248527 0.562    
7   4.57e-18  0.2471 8.526 0.007855 0.563   *
8  -8.97e-01 -1.9019 0.287 0.319219 0.562    
9  -1.91e-01  0.4042 7.579 0.020741 0.562   *
10 -4.77e-17 -1.2890 1.582 0.177221 0.562    
11  5.24e-17 -0.4042 7.579 0.020741 0.562   *
12  6.08e-01  1.2890 1.582 0.177221 0.562    
13  1.31e-01 -0.2783 8.366 0.009941 0.562   *
14  8.13e-17 -2.1970 0.115 0.384880 0.563   *
15  3.60e-17  0.2783 8.366 0.009941 0.562   *
16  1.04e+00  2.1970 0.115 0.384880 0.562   *
influenceIndexPlot(modelo.dcl)

Cumplimiento de supuestos del modelo lineal general

Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dcl,
                 reps = 5000,
                 max.lag = 5)
 lag Autocorrelation D-W Statistic p-value
   1     -0.06276847      1.956937  0.6012
   2     -0.59793814      2.894008  0.0600
   3      0.10266323      1.323561  0.6764
   4      0.24774485      0.929768  0.0164
   5     -0.08574957      1.519008  0.4408
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dcl, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dcl
DW = 1.9569, p-value = 0.6117
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dcl))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dcl)
W = 0.95851, p-value = 0.6348
ad.test(rstudent(modelo.dcl))

    Anderson-Darling normality test

data:  rstudent(modelo.dcl)
A = 0.29833, p-value = 0.5438
lillie.test(rstudent(modelo.dcl))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dcl)
D = 0.12781, p-value = 0.6915
ks.test(rstudent(modelo.dcl), "pnorm",
        alternative = "two.sided")

    Exact one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dcl)
D = 0.12219, p-value = 0.947
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dcl))

    Cramer-von Mises normality test

data:  rstudent(modelo.dcl)
W = 0.051311, p-value = 0.4692
pearson.test(rstudent(modelo.dcl))

    Pearson chi-square normality test

data:  rstudent(modelo.dcl)
P = 2.375, p-value = 0.6671
sf.test(rstudent(modelo.dcl))

    Shapiro-Francia normality test

data:  rstudent(modelo.dcl)
W = 0.9681, p-value = 0.7196

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento

ncvTest(modelo.dcl)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.2227318, Df = 1, p = 0.63697
bptest(modelo.dcl)

    studentized Breusch-Pagan test

data:  modelo.dcl
BP = 14.686, df = 8, p-value = 0.06555
bptest(modelo.dcl, studentize = F)

    Breusch-Pagan test

data:  modelo.dcl
BP = 7.5483, df = 8, p-value = 0.4788
olsrr::ols_test_breusch_pagan(modelo.dcl)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

              Data               
 --------------------------------
 Response : rdto 
 Variables: fitted values of rdto 

        Test Summary         
 ----------------------------
 DF            =    1 
 Chi2          =    0.2227318 
 Prob > Chi2   =    0.6369663

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que no se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dcl %>% gvlma()

Call:
lm(formula = rdto ~ sp + trt + row + col, data = data)

Coefficients:
(Intercept)        spsp2         trtB         row2         row3         row4  
      8.050        0.900       -3.625        2.100        1.350        0.800  
       col2         col3         col4  
      0.650       -0.350        0.750  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                       Value p-value                Decision
Global Stat        1.180e+00  0.8814 Assumptions acceptable.
Skewness           5.756e-05  0.9939 Assumptions acceptable.
Kurtosis           6.299e-01  0.4274 Assumptions acceptable.
Link Function      1.635e-01  0.6860 Assumptions acceptable.
Heteroscedasticity 3.864e-01  0.5342 Assumptions acceptable.

Análisis de varianza

\[Y_{ijkl} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_{l} + \epsilon_{ijkl}\]

\[\hat{Y}_{ijkl} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_{l}\] Dónde:

\(Y_{ijkl}\) = Valor observado de la variable respuesta.

\(\hat{Y}_{ijkl}\) = Valor ajustado de la variable respuesta.

\(\mu\) = Promedio observado de la variable respuesta.

\(\tau_{i}\) = Efecto del i-ésimo nivel del factor A.

\(\beta_{j}\) = Efecto del j-ésimo nivel del factor B.

\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor Fila.

\(\delta_{l}\) = Efecto del k-ésimo nivel del factor Columna.

\(\epsilon_{ijkl}\) = Residuo observado del modelo.

Pruebas de hipótesis

Para el factor A (Especie):

\(H_0: \tau_{A1} = \tau_{A2} = 0\)

\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)

Para el factor B (Tratamiento):

\(H_0: \beta_{B1} = \beta_{B2} = 0\)

\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_{j} \neq 0\text{; en al menos un nivel del factor B.}\)

anova(modelo.dcl, test = "F")
Analysis of Variance Table

Response: rdto
          Df Sum Sq Mean Sq F value    Pr(>F)    
sp         1  3.240   3.240  1.9485 0.2054287    
trt        1 52.563  52.563 31.6098 0.0007972 ***
row        3  9.428   3.143  1.8898 0.2196546    
col        3  3.327   1.109  0.6670 0.5985006    
Residuals  7 11.640   1.663                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 1, 6)
[1] 5.987378

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 1, 6)
[1] 5.987378

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor A tienen un efecto sobre el rendimiento estadísticamente similar a 0.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor B.

agricolae::cv.model(modelo.dcl)
[1] 16.09383

Comparaciones de medias

Para los niveles del factor B:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

Prueba de LSD

Para los niveles del Factor B:

agricolae::LSD.test(modelo.dcl, trt = "trt", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

LSD t Test for rdto 

Mean Square Error:  1.662857 

trt,  means and individual ( 95 %) CI

   rdto      std r      LCL       UCL Min  Max
A 9.825 1.404839 8 8.746936 10.903064   8 12.0
B 6.200 1.405093 8 5.121936  7.278064   5  9.2

Alpha: 0.05 ; DF Error: 7
Critical Value of t: 2.364624 

least Significant Difference: 1.524613 

Treatments with the same letter are not significantly different.

   rdto groups
A 9.825      a
B 6.200      b
agricolae::LSD.test(modelo.dcl, trt = "trt", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "trt"

LSD t Test for rdto 

Mean Square Error:  1.662857 

trt,  means and individual ( 95 %) CI

   rdto      std r      LCL       UCL Min  Max
A 9.825 1.404839 8 8.746936 10.903064   8 12.0
B 6.200 1.405093 8 5.121936  7.278064   5  9.2

Alpha: 0.05 ; DF Error: 7
Critical Value of t: 2.364624 

Comparison between treatments means

      difference pvalue signif.      LCL      UCL
A - B      3.625  8e-04     *** 2.100387 5.149613

Diseño cuadrado grecolatino

Planeamiento

Crear un libro de campo con el paquete agricolae

T1<-c("a","b","c","d")
T2<-c("v","w","x","y")

salida <- agricolae::design.graeco(trt1 = T1,
                                   trt2 = T2,
                                   serie = 3,
                                   seed = 123,
                                   kinds = "Super-Duper",
                                   randomization = TRUE)
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
salida$sketch
     [,1]  [,2]  [,3]  [,4] 
[1,] "c v" "a x" "d w" "b y"
[2,] "a w" "c y" "b v" "d x"
[3,] "d y" "b w" "c x" "a v"
[4,] "b x" "d v" "a y" "c w"
print(matrix(as.numeric(salida$book[,1]),
             byrow=TRUE,
             ncol=4))
     [,1] [,2] [,3] [,4]
[1,] 1001 1002 1003 1004
[2,] 2001 2002 2003 2004
[3,] 3001 3002 3003 3004
[4,] 4001 4002 4003 4004

Guardar el libro generado

write.table(salida$book,
            "books/graeco.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida$book,
           "books/graeco.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)

Guardar el libro de campo generado

write.table(salida %>% zigzag(),
            "books/graeco.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida %>% zigzag(),
           "books/graeco.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
agricolaeplotr::plot_graeco(
  design = salida,
  factor_name = "T1",
  reverse_y = TRUE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas")
agricolaeplotr::plot_graeco(
  design = salida,
  factor_name = "T2",
  reverse_y = TRUE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Columnas",
       y = "Filas")

Crear un libro de campo con el paquete edibble

menu_graeco()
design("Graeco-Latin Square Design") %>%
  set_units(row = 8,
            col = 8,
            unit = crossed_by(row, col)) %>%
  set_trts(trt1 = 8,
           trt2 = 8) %>%
  allot_trts(trt1 ~ unit,
             trt2 ~ unit) %>%
  assign_trts("random", seed = 802) %>%
  serve_table()
graeco <- takeout(menu_graeco(t = c(4,4),
                              seed = 123))
graeco %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
graeco2 <- design("Graeco-Latin Square Design") %>%
  set_units(row = 16,
            col = 16,
            unit = crossed_by(row, col)) %>%
  set_trts(trt1 = c("Amazon", "Coari x Lame", "Coari x Yangambí", "Criollo"),
           trt2 = c(0, 30, 60, 90)) %>%
  allot_trts(trt1 + trt2 ~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
graeco2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
write.table(graeco2 %>% as.data.frame(),
            "books/graeco.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(graeco2 %>% as.data.frame(),
           "books/graeco.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
deggust::autoplot(lsd2)
plot(graeco2)

Análisis de DCGL

Importación de datos

archivos <- list.files(pattern = "DCGL.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor) %>%
  mutate(row = factor(`F`),
         col = factor(C)) %>%
  dplyr::select(-c(`F`,C))

Creación del modelo lineal con interacción

modelo.dcgl <- lm(rdto ~ Variedad + Densidad + row + col, data = data)

Definición del modelo

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*A_3 + \beta_3*A_4 + \beta_4*A_5 + \beta_5*B_2 + \beta_6*B_3 + \beta_7*B_4 + \beta_8*B_5 + \beta_9*F_{2} + \beta_{10}*F_{3} + \beta_{11}*F_{4} + \beta_{12}*F_{5} + \beta_{13}*C_2 + \beta_{14}*C_3 + \beta_{15}*C_4 + \beta_{16}*C_5 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*A_3 + \beta_3*A_4 + \beta_4*A_5 + \beta_5*B_2 + \beta_6*B_3 + \beta_7*B_4 + \beta_8*B_5 + \beta_9*F_{2} + \beta_{10}*F_{3} + \beta_{11}*F_{4} + \beta_{12}*F_{5} + \beta_{13}*C_2 + \beta_{14}*C_3 + \beta_{15}*C_4 + \beta_{16}*C_5\]

summary(modelo.dcgl)

Call:
lm(formula = rdto ~ Variedad + Densidad + row + col, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
  -2.8   -1.0    0.2    1.0    2.0 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   24.400      1.994  12.234 1.85e-06 ***
VariedadB     -8.000      1.530  -5.230 0.000793 ***
VariedadC     -4.800      1.530  -3.138 0.013850 *  
VariedadD     -8.600      1.530  -5.622 0.000497 ***
VariedadE    -10.600      1.530  -6.929 0.000121 ***
Densidadb      0.400      1.530   0.261 0.800322    
Densidadc      1.600      1.530   1.046 0.326154    
Densidadd     -0.200      1.530  -0.131 0.899206    
Densidade      1.200      1.530   0.784 0.455366    
row2          -0.200      1.530  -0.131 0.899206    
row3          -0.800      1.530  -0.523 0.615162    
row4          -1.400      1.530  -0.915 0.386837    
row5          -1.600      1.530  -1.046 0.326154    
col2          -0.200      1.530  -0.131 0.899206    
col3           0.600      1.530   0.392 0.705129    
col4          -1.200      1.530  -0.784 0.455366    
col5          -2.200      1.530  -1.438 0.188326    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.419 on 8 degrees of freedom
Multiple R-squared:  0.8927,    Adjusted R-squared:  0.678 
F-statistic: 4.158 on 16 and 8 DF,  p-value: 0.02355

Verificación visual de los supuestos del modelo

performance::check_model(modelo.dcgl)
ggResidpanel::resid_panel(modelo.dcgl)
influence.measures(modelo.dcgl)
Influence measures of
     lm(formula = rdto ~ Variedad + Densidad + row + col, data = data) :

    dfb.1_  dfb.VrdB  dfb.VrdC  dfb.VrdD  dfb.VrdE dfb.Dnsddb dfb.Dnsddc
1   1.7513 -6.72e-01 -6.72e-01 -6.72e-01 -6.72e-01  -6.72e-01  -6.72e-01
2   0.0235  7.65e-02  3.15e-18 -1.74e-17 -7.73e-18  -1.06e-17   7.65e-02
3   0.0000  0.00e+00  0.00e+00  0.00e+00  0.00e+00   0.00e+00   0.00e+00
4   0.0235 -2.66e-18  0.00e+00  7.65e-02  1.81e-17   7.65e-02  -2.74e-17
5  -0.2739 -3.02e-17  1.16e-17  6.11e-16 -8.93e-01   1.38e-16  -1.92e-16
6  -0.0712 -2.32e-01  2.67e-17 -8.12e-17  2.55e-17  -2.32e-01   8.63e-17
7  -0.4109  4.34e-16  8.93e-01  2.39e-16  4.70e-17   1.27e-16   3.84e-16
8  -0.4756  7.31e-16  5.55e-16 -1.55e+00  2.18e-16   1.55e+00   1.55e+00
9  -0.2220  1.55e-16 -8.14e-17 -3.67e-17  4.82e-01  -3.35e-17   4.82e-01
10  0.0235 -7.65e-02 -7.65e-02 -7.65e-02 -7.65e-02  -1.91e-17  -2.62e-17
11 -0.4756  5.38e-17 -1.55e+00 -7.74e-16  1.90e-16  -1.08e-16  -1.55e+00
12 -0.1068  6.83e-17 -1.66e-17  2.32e-01  1.43e-17  -2.84e-17   7.34e-17
13 -0.4109 -7.45e-17 -2.65e-17 -2.71e-16  8.93e-01   8.93e-01   5.26e-17
14 -0.1761  5.74e-01  5.74e-01  5.74e-01  5.74e-01   4.05e-17   8.23e-17
15  0.2060  6.72e-01 -2.81e-17 -9.32e-17  5.60e-17  -6.72e-01  -6.72e-01
16  0.2060 -2.47e-16 -2.41e-16  6.72e-01 -2.30e-16  -1.20e-17  -4.81e-17
17 -0.1761  7.98e-17  1.52e-16  8.73e-17 -5.74e-01   5.74e-01   5.74e-01
18  0.1214 -3.96e-01 -3.96e-01 -3.96e-01 -3.96e-01   9.61e-17   3.96e-01
19  0.1821 -3.96e-01  5.04e-17  4.49e-17  1.91e-17  -1.34e-17  -4.84e-17
20  0.0352  1.33e-17 -7.65e-02  2.45e-17  1.14e-17  -7.65e-02  -1.20e-17
21  0.0235 -3.36e-17 -2.90e-17  1.16e-17  7.65e-02  -3.25e-18  -1.68e-17
22 -0.1761  5.74e-01  5.74e-01  5.74e-01  5.74e-01  -5.74e-01   3.20e-17
23  0.0352 -7.65e-02 -1.59e-17 -2.03e-17  5.04e-18  -1.75e-17  -2.74e-17
24  0.1214 -2.88e-16  3.96e-01 -1.64e-16 -7.99e-17  -3.96e-01  -3.96e-01
25 -0.0707  6.67e-17  1.10e-17  1.54e-01  2.83e-17   0.00e+00   1.54e-01
   dfb.Dnsddd dfb.Densdd  dfb.row2  dfb.row3  dfb.row4  dfb.row5  dfb.col2
1   -6.72e-01  -6.72e-01 -6.72e-01 -6.72e-01 -6.72e-01 -6.72e-01 -6.72e-01
2    2.08e-17   3.69e-18  7.65e-02  1.54e-17  3.66e-18 -3.22e-17 -7.65e-02
3    0.00e+00   0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00
4   -1.83e-17  -1.68e-18  8.14e-18  2.19e-17  7.65e-02 -2.02e-18 -7.65e-02
5   -8.93e-01  -2.86e-16  1.17e-16  2.71e-16  8.14e-16 -8.93e-01  8.93e-01
6    1.70e-17   8.96e-17  2.32e-01  2.32e-01  2.32e-01  2.32e-01 -2.32e-01
7    8.93e-01   6.27e-17  8.93e-01  8.78e-17  4.98e-16  5.96e-16  8.93e-01
8    1.55e+00   1.55e+00  2.15e-16 -1.55e+00  6.47e-16 -1.36e-17 -1.55e+00
9   -9.98e-17  -5.29e-17 -9.34e-17 -1.73e-16  4.82e-01 -1.23e-16  4.82e-01
10  -2.91e-17   7.65e-02 -1.37e-17  5.15e-18 -1.75e-17  7.65e-02  7.65e-02
11  -2.16e-17   3.81e-16  1.55e+00  1.55e+00  1.55e+00  1.55e+00 -5.38e-17
12   7.06e-17   2.32e-01  2.32e-01 -2.18e-17 -2.36e-18 -4.18e-17  6.44e-17
13   1.08e-16   1.21e-16  2.48e-16  8.93e-01  1.31e-16  2.90e-16  3.10e-17
14  -5.74e-01   3.27e-17  1.70e-16  4.11e-17 -5.74e-01  1.51e-17  2.59e-16
15  -6.72e-01  -6.72e-01  2.21e-16  2.12e-16  2.04e-16  6.72e-01 -2.33e-17
16   6.72e-01  -1.50e-16 -6.72e-01 -6.72e-01 -6.72e-01 -6.72e-01  3.63e-16
17   5.74e-01   5.74e-01 -5.74e-01 -2.50e-17 -1.10e-16 -1.16e-16  9.96e-17
18   1.31e-16   8.85e-17  1.37e-16  3.96e-01  1.71e-16  1.49e-16 -2.84e-16
19  -6.01e-17  -3.96e-01 -2.13e-17 -8.51e-17 -3.96e-01 -7.99e-17 -2.01e-16
20  -1.38e-17  -2.18e-17  2.04e-17  2.07e-17  2.91e-17 -7.65e-02 -2.14e-17
21  -1.75e-17   7.65e-02 -7.65e-02 -7.65e-02 -7.65e-02 -7.65e-02 -7.28e-18
22  -4.17e-17  -1.18e-16 -5.74e-01  6.67e-17  1.06e-16  1.16e-16 -1.04e-16
23  -7.65e-02  -1.85e-17 -2.12e-17 -7.65e-02 -2.43e-17 -4.70e-18  6.29e-18
24  -3.96e-01  -3.96e-01 -1.10e-16 -5.67e-17  3.96e-01 -6.25e-17  1.20e-16
25   4.63e-17   3.31e-17  1.11e-17 -1.56e-17 -4.67e-17  1.54e-01  9.84e-17
    dfb.col3  dfb.col4  dfb.col5     dffit    cov.r   cook.d  hat inf
1  -6.72e-01 -6.72e-01 -6.72e-01  1.75e+00 1.25e+00 1.71e-01 0.68   *
2  -7.65e-02 -7.65e-02 -7.65e-02  2.00e-01 2.89e+01 2.67e-03 0.68   *
3   0.00e+00  0.00e+00  0.00e+00  1.16e-15 3.02e+01 9.07e-32 0.68   *
4  -7.65e-02 -7.65e-02 -7.65e-02  2.00e-01 2.89e+01 2.67e-03 0.68   *
5   8.93e-01  8.93e-01  8.93e-01 -2.33e+00 1.54e-01 2.67e-01 0.68    
6  -5.64e-17 -2.58e-17 -5.60e-17 -6.05e-01 2.00e+01 2.40e-02 0.68   *
7  -5.58e-18 -7.11e-17  7.84e-17  2.33e+00 1.54e-01 2.67e-01 0.68    
8   5.96e-17  4.02e-16  6.80e-17 -4.04e+00 1.02e-04 5.24e-01 0.68   *
9   1.91e-16  8.22e-17  6.35e-17  1.26e+00 5.42e+00 9.62e-02 0.68    
10 -2.15e-17 -2.73e-17 -2.69e-17  2.00e-01 2.89e+01 2.67e-03 0.68   *
11 -1.55e+00  1.91e-16  0.00e+00 -4.04e+00 1.02e-04 5.24e-01 0.68   *
12  2.32e-01  2.33e-17  5.09e-17  6.05e-01 2.00e+01 2.40e-02 0.68   *
13  8.93e-01  2.28e-16  2.35e-16  2.33e+00 1.54e-01 2.67e-01 0.68    
14 -5.74e-01  1.17e-16  2.01e-16 -1.50e+00 2.79e+00 1.31e-01 0.68    
15  6.72e-01 -3.93e-17 -1.47e-17  1.75e+00 1.25e+00 1.71e-01 0.68    
16  2.41e-16  6.72e-01  2.65e-16  1.75e+00 1.25e+00 1.71e-01 0.68    
17  0.00e+00 -5.74e-01  1.01e-16 -1.50e+00 2.79e+00 1.31e-01 0.68    
18 -3.69e-16  3.96e-01 -2.43e-16  1.03e+00 9.34e+00 6.68e-02 0.68   *
19 -2.55e-16 -3.96e-01 -1.30e-16 -1.03e+00 9.34e+00 6.68e-02 0.68   *
20 -2.19e-17 -7.65e-02 -2.02e-17 -2.00e-01 2.89e+01 2.67e-03 0.68   *
21 -2.23e-17 -1.75e-17  7.65e-02  2.00e-01 2.89e+01 2.67e-03 0.68   *
22 -5.57e-17 -4.36e-17 -5.74e-01 -1.50e+00 2.79e+00 1.31e-01 0.68    
23  3.10e-17 -5.82e-18 -7.65e-02 -2.00e-01 2.89e+01 2.67e-03 0.68   *
24  1.64e-16  1.50e-16  3.96e-01  1.03e+00 9.34e+00 6.68e-02 0.68   *
25  6.74e-17  7.01e-17  1.54e-01  4.01e-01 2.52e+01 1.07e-02 0.68   *
influenceIndexPlot(modelo.dcgl)

Cumplimiento de supuestos del modelo lineal general

Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dcgl,
                 reps = 5000,
                 max.lag = 5)
 lag Autocorrelation D-W Statistic p-value
   1      -0.3547009     2.6512821  0.3528
   2      -0.2170940     2.3538462  0.9264
   3       0.5589744     0.8008547  0.0700
   4      -0.4871795     2.8504274  0.0612
   5      -0.1290598     2.0478632  0.9248
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dcgl, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dcgl
DW = 2.6513, p-value = 0.3357
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dcgl))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dcgl)
W = 0.92827, p-value = 0.07926
ad.test(rstudent(modelo.dcgl))

    Anderson-Darling normality test

data:  rstudent(modelo.dcgl)
A = 0.53976, p-value = 0.1497
lillie.test(rstudent(modelo.dcgl))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dcgl)
D = 0.15131, p-value = 0.1473
ks.test(rstudent(modelo.dcgl), "pnorm",
        alternative = "two.sided")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dcgl)
D = 0.12555, p-value = 0.8256
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dcgl))

    Cramer-von Mises normality test

data:  rstudent(modelo.dcgl)
W = 0.07781, p-value = 0.2112
pearson.test(rstudent(modelo.dcgl))

    Pearson chi-square normality test

data:  rstudent(modelo.dcgl)
P = 4.76, p-value = 0.4459
sf.test(rstudent(modelo.dcgl))

    Shapiro-Francia normality test

data:  rstudent(modelo.dcgl)
W = 0.93314, p-value = 0.09472

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento

ncvTest(modelo.dcgl)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.02393361, Df = 1, p = 0.87705
bptest(modelo.dcgl)

    studentized Breusch-Pagan test

data:  modelo.dcgl
BP = 14.799, df = 16, p-value = 0.5394
bptest(modelo.dcgl, studentize = F)

    Breusch-Pagan test

data:  modelo.dcgl
BP = 10.125, df = 16, p-value = 0.86
olsrr::ols_test_breusch_pagan(modelo.dcgl)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

              Data               
 --------------------------------
 Response : rdto 
 Variables: fitted values of rdto 

        Test Summary          
 -----------------------------
 DF            =    1 
 Chi2          =    0.02393361 
 Prob > Chi2   =    0.877054

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que no se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dcgl %>% gvlma()

Call:
lm(formula = rdto ~ Variedad + Densidad + row + col, data = data)

Coefficients:
(Intercept)    VariedadB    VariedadC    VariedadD    VariedadE    Densidadb  
       24.4         -8.0         -4.8         -8.6        -10.6          0.4  
  Densidadc    Densidadd    Densidade         row2         row3         row4  
        1.6         -0.2          1.2         -0.2         -0.8         -1.4  
       row5         col2         col3         col4         col5  
       -1.6         -0.2          0.6         -1.2         -2.2  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                    Value p-value                   Decision
Global Stat        6.0111 0.19832    Assumptions acceptable.
Skewness           0.8570 0.35458    Assumptions acceptable.
Kurtosis           0.4157 0.51910    Assumptions acceptable.
Link Function      3.9267 0.04752 Assumptions NOT satisfied!
Heteroscedasticity 0.8117 0.36761    Assumptions acceptable.

Análisis de varianza

\[Y_{ijkl} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_{l} + \epsilon_{ijkl}\]

\[\hat{Y}_{ijkl} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_{l}\] Dónde:

\(Y_{ijkl}\) = Valor observado de la variable respuesta.

\(\hat{Y}_{ijkl}\) = Valor ajustado de la variable respuesta.

\(\mu\) = Promedio observado de la variable respuesta.

\(\tau_{i}\) = Efecto del i-ésimo nivel del factor A.

\(\beta_{j}\) = Efecto del j-ésimo nivel del factor B.

\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor Fila.

\(\delta_{l}\) = Efecto del l-ésimo nivel del factor Columna.

\(\epsilon_{ijkl}\) = Residuo observado del modelo.

Pruebas de hipótesis

Para el factor A (Especie):

\(H_0: \tau_{A1} = \tau_{A2} = 0\)

\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)

Para el factor B (Tratamiento):

\(H_0: \beta_{B1} = \beta_{B2} = 0\)

\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_{j} \neq 0\text{; en al menos un nivel del factor B.}\)

anova(modelo.dcgl, test = "F")
Analysis of Variance Table

Response: rdto
          Df Sum Sq Mean Sq F value   Pr(>F)    
Variedad   4  342.8   85.70 14.6496 0.000941 ***
Densidad   4   12.0    3.00  0.5128 0.728900    
row        4   10.0    2.50  0.4274 0.785447    
col        4   24.4    6.10  1.0427 0.442543    
Residuals  8   46.8    5.85                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 4, 8)
[1] 3.837853

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 4, 8)
[1] 3.837853

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, l menos un nivel del factor A tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor A.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor B tienen un efecto sobre el rendimiento estadísticamente similar a 0.

agricolae::cv.model(modelo.dcgl)
[1] 14.06208

Comparaciones de medias

Para los niveles del factor A:

  • A1 vs A2:

\(H_0: \mu_{A1} - \mu_{A2} = 0\)

\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)

  • A1 vs A3:

\(H_0: \mu_{A1} - \mu_{A3} = 0\)

\(H_1: \mu_{A1} - \mu_{A3} \neq 0\)

  • A1 vs A4:

\(H_0: \mu_{A1} - \mu_{A4} = 0\)

\(H_1: \mu_{A1} - \mu_{A4} \neq 0\)

  • A1 vs A5:

\(H_0: \mu_{A1} - \mu_{A5} = 0\)

\(H_1: \mu_{A1} - \mu_{A5} \neq 0\)

  • A2 vs A3:

\(H_0: \mu_{A2} - \mu_{A3} = 0\)

\(H_1: \mu_{A2} - \mu_{A3} \neq 0\)

  • A2 vs A4:

\(H_0: \mu_{A2} - \mu_{A4} = 0\)

\(H_1: \mu_{A2} - \mu_{A4} \neq 0\)

  • A2 vs A5:

\(H_0: \mu_{A2} - \mu_{A5} = 0\)

\(H_1: \mu_{A2} - \mu_{A5} \neq 0\)

  • A3 vs A4:

\(H_0: \mu_{A3} - \mu_{A4} = 0\)

\(H_1: \mu_{A3} - \mu_{A4} \neq 0\)

  • A3 vs A5:

\(H_0: \mu_{A3} - \mu_{A5} = 0\)

\(H_1: \mu_{A3} - \mu_{A5} \neq 0\)

  • A4 vs A5:

\(H_0: \mu_{A4} - \mu_{A5} = 0\)

\(H_1: \mu_{A4} - \mu_{A5} \neq 0\)

Prueba de LSD

Para los niveles del Factor A:

agricolae::LSD.test(modelo.dcgl, trt = "Variedad", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dcgl ~ "Variedad"

LSD t Test for rdto 

Mean Square Error:  5.85 

Variedad,  means and individual ( 95 %) CI

  rdto      std r      LCL      UCL Min Max
A 23.6 2.073644 5 21.10568 26.09432  21  26
B 15.6 2.073644 5 13.10568 18.09432  13  18
C 18.8 1.788854 5 16.30568 21.29432  17  21
D 15.0 2.236068 5 12.50568 17.49432  12  18
E 13.0 2.549510 5 10.50568 15.49432  10  16

Alpha: 0.05 ; DF Error: 8
Critical Value of t: 2.306004 

least Significant Difference: 3.527508 

Treatments with the same letter are not significantly different.

  rdto groups
A 23.6      a
C 18.8      b
B 15.6     bc
D 15.0      c
E 13.0      c
agricolae::LSD.test(modelo.dcgl, trt = "Variedad", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dcgl ~ "Variedad"

LSD t Test for rdto 

Mean Square Error:  5.85 

Variedad,  means and individual ( 95 %) CI

  rdto      std r      LCL      UCL Min Max
A 23.6 2.073644 5 21.10568 26.09432  21  26
B 15.6 2.073644 5 13.10568 18.09432  13  18
C 18.8 1.788854 5 16.30568 21.29432  17  21
D 15.0 2.236068 5 12.50568 17.49432  12  18
E 13.0 2.549510 5 10.50568 15.49432  10  16

Alpha: 0.05 ; DF Error: 8
Critical Value of t: 2.306004 

Comparison between treatments means

      difference pvalue signif.       LCL       UCL
A - B        8.0 0.0008     ***  4.472492 11.527508
A - C        4.8 0.0138       *  1.272492  8.327508
A - D        8.6 0.0005     ***  5.072492 12.127508
A - E       10.6 0.0001     ***  7.072492 14.127508
B - C       -3.2 0.0698       . -6.727508  0.327508
B - D        0.6 0.7051         -2.927508  4.127508
B - E        2.6 0.1276         -0.927508  6.127508
C - D        3.8 0.0379       *  0.272492  7.327508
C - E        5.8 0.0053      **  2.272492  9.327508
D - E        2.0 0.2274         -1.527508  5.527508

Transformación de datos

Transformación logarítmica

data <- data %>%
  mutate(rdto_logy_1 = log(rdto #+ 1
                           ),
         rdto_logiy = log(1/rdto))
modelo.dcgl_log <- lm(rdto_logy_1 ~ Variedad + Densidad + row + col, data = data)
summary(modelo.dcgl_log)

Call:
lm(formula = rdto_logy_1 ~ Variedad + Densidad + row + col, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.18885 -0.07671  0.01741  0.06225  0.12955 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.167131   0.130803  24.213 9.03e-09 ***
VariedadB   -0.418082   0.100321  -4.167 0.003133 ** 
VariedadC   -0.227894   0.100321  -2.272 0.052755 .  
VariedadD   -0.459123   0.100321  -4.577 0.001810 ** 
VariedadE   -0.608904   0.100321  -6.070 0.000299 ***
Densidadb    0.064728   0.100321   0.645 0.536856    
Densidadc    0.118694   0.100321   1.183 0.270720    
Densidadd   -0.006327   0.100321  -0.063 0.951257    
Densidade    0.092168   0.100321   0.919 0.385101    
row2        -0.008996   0.100321  -0.090 0.930754    
row3        -0.055122   0.100321  -0.549 0.597693    
row4        -0.066650   0.100321  -0.664 0.525140    
row5        -0.105308   0.100321  -1.050 0.324526    
col2         0.006657   0.100321   0.066 0.948721    
col3         0.075364   0.100321   0.751 0.474033    
col4        -0.058641   0.100321  -0.585 0.574972    
col5        -0.101629   0.100321  -1.013 0.340699    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1586 on 8 degrees of freedom
Multiple R-squared:  0.8654,    Adjusted R-squared:  0.5962 
F-statistic: 3.214 on 16 and 8 DF,  p-value: 0.04945
performance::check_model(modelo.dcgl_log)
ggResidpanel::resid_panel(modelo.dcgl_log)

Estadísticas globales

modelo.dcgl_log %>% gvlma()

Call:
lm(formula = rdto_logy_1 ~ Variedad + Densidad + row + col, data = data)

Coefficients:
(Intercept)    VariedadB    VariedadC    VariedadD    VariedadE    Densidadb  
   3.167131    -0.418082    -0.227894    -0.459123    -0.608904     0.064728  
  Densidadc    Densidadd    Densidade         row2         row3         row4  
   0.118694    -0.006327     0.092168    -0.008996    -0.055122    -0.066650  
       row5         col2         col3         col4         col5  
  -0.105308     0.006657     0.075364    -0.058641    -0.101629  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                     Value   p-value                   Decision
Global Stat        19.7030 5.715e-04 Assumptions NOT satisfied!
Skewness            0.8430 3.585e-01    Assumptions acceptable.
Kurtosis            0.3839 5.355e-01    Assumptions acceptable.
Link Function      17.6976 2.589e-05 Assumptions NOT satisfied!
Heteroscedasticity  0.7785 3.776e-01    Assumptions acceptable.

Transformación de Box y Cox

set.seed(123)
p <- car::boxCox(object = modelo.dcgl, lambda = seq(-2,2, 0.05))
lambda <- p$x[which(p$y == max(p$y))]
lambda
[1] 1.474747
data <- data %>%
  mutate(rdto_bc = ((rdto^lambda)-1)/lambda)
modelo.dcgl_bc <- lm(rdto_bc ~ Variedad + Densidad + row + col, data = data)
summary(modelo.dcgl_bc)

Call:
lm(formula = rdto_bc ~ Variedad + Densidad + row + col, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.7052  -3.3492   0.6157   4.1061   8.0309 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  75.6141     7.5012  10.080 8.00e-06 ***
VariedadB   -32.7583     5.7531  -5.694 0.000458 ***
VariedadC   -20.4750     5.7531  -3.559 0.007412 ** 
VariedadD   -34.9132     5.7531  -6.069 0.000300 ***
VariedadE   -41.8203     5.7531  -7.269 8.64e-05 ***
Densidadb     0.2679     5.7531   0.047 0.964005    
Densidadc     5.4774     5.7531   0.952 0.368937    
Densidadd    -1.0843     5.7531  -0.188 0.855202    
Densidade     3.9837     5.7531   0.692 0.508257    
row2         -1.0465     5.7531  -0.182 0.860191    
row3         -2.9093     5.7531  -0.506 0.626716    
row4         -5.9747     5.7531  -1.039 0.329401    
row5         -6.0225     5.7531  -1.047 0.325780    
col2         -1.2916     5.7531  -0.225 0.827993    
col3          1.1026     5.7531   0.192 0.852794    
col4         -4.9637     5.7531  -0.863 0.413376    
col5         -9.3407     5.7531  -1.624 0.143118    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.096 on 8 degrees of freedom
Multiple R-squared:  0.9019,    Adjusted R-squared:  0.7058 
F-statistic: 4.598 on 16 and 8 DF,  p-value: 0.01735
performance::check_model(modelo.dcgl_bc)
ggResidpanel::resid_panel(modelo.dcgl_bc)

Estadísticas globales

modelo.dcgl_bc %>% gvlma()

Call:
lm(formula = rdto_bc ~ Variedad + Densidad + row + col, data = data)

Coefficients:
(Intercept)    VariedadB    VariedadC    VariedadD    VariedadE    Densidadb  
    75.6141     -32.7583     -20.4750     -34.9132     -41.8203       0.2679  
  Densidadc    Densidadd    Densidade         row2         row3         row4  
     5.4774      -1.0843       3.9837      -1.0465      -2.9093      -5.9747  
       row5         col2         col3         col4         col5  
    -6.0225      -1.2916       1.1026      -4.9637      -9.3407  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                      Value p-value                Decision
Global Stat        2.073858  0.7222 Assumptions acceptable.
Skewness           0.770212  0.3802 Assumptions acceptable.
Kurtosis           0.416263  0.5188 Assumptions acceptable.
Link Function      0.002935  0.9568 Assumptions acceptable.
Heteroscedasticity 0.884448  0.3470 Assumptions acceptable.

Análisis de varianza

anova(modelo.dcgl_bc, test = "F")
Analysis of Variance Table

Response: rdto_bc
          Df Sum Sq Mean Sq F value    Pr(>F)    
Variedad   4 5409.6 1352.41 16.3441 0.0006449 ***
Densidad   4  160.9   40.22  0.4860 0.7464677    
row        4  153.1   38.28  0.4627 0.7619143    
col        4  363.8   90.95  1.0991 0.4194963    
Residuals  8  662.0   82.75                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 4, 8)
[1] 3.837853

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 4, 8)
[1] 3.837853

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, l menos un nivel del factor A tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor A.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor B tienen un efecto sobre el rendimiento estadísticamente similar a 0.

agricolae::cv.model(modelo.dcgl_bc)
[1] 20.09809

Comparaciones de medias

Prueba de LSD

Para los niveles del Factor A:

agricolae::LSD.test(modelo.dcgl_bc, trt = "Variedad", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dcgl_bc ~ "Variedad"

LSD t Test for rdto_bc 

Mean Square Error:  82.74609 

Variedad,  means and individual ( 95 %) CI

   rdto_bc      std r      LCL      UCL      Min      Max
A 71.25381 9.270546 5 61.87281 80.63480 59.74765 82.11818
B 38.49554 7.601826 5 29.11455 47.87653 29.11171 47.46039
C 50.77879 7.205075 5 41.39780 60.15978 43.56891 59.74765
D 36.34060 8.079757 5 26.95961 45.72160 25.79486 47.46039
E 29.43354 8.603178 5 20.05255 38.81454 19.55350 39.78465

Alpha: 0.05 ; DF Error: 8
Critical Value of t: 2.306004 

least Significant Difference: 13.26673 

Treatments with the same letter are not significantly different.

   rdto_bc groups
A 71.25381      a
C 50.77879      b
B 38.49554     bc
D 36.34060      c
E 29.43354      c
agricolae::LSD.test(modelo.dcgl_bc, trt = "Variedad", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dcgl_bc ~ "Variedad"

LSD t Test for rdto_bc 

Mean Square Error:  82.74609 

Variedad,  means and individual ( 95 %) CI

   rdto_bc      std r      LCL      UCL      Min      Max
A 71.25381 9.270546 5 61.87281 80.63480 59.74765 82.11818
B 38.49554 7.601826 5 29.11455 47.87653 29.11171 47.46039
C 50.77879 7.205075 5 41.39780 60.15978 43.56891 59.74765
D 36.34060 8.079757 5 26.95961 45.72160 25.79486 47.46039
E 29.43354 8.603178 5 20.05255 38.81454 19.55350 39.78465

Alpha: 0.05 ; DF Error: 8
Critical Value of t: 2.306004 

Comparison between treatments means

      difference pvalue signif.        LCL        UCL
A - B  32.758269 0.0005     ***  19.491541 46.0249967
A - C  20.475019 0.0074      **   7.208291 33.7417469
A - D  34.913205 0.0003     ***  21.646478 48.1799331
A - E  41.820263 0.0001     ***  28.553535 55.0869909
B - C -12.283250 0.0653       . -25.549978  0.9834778
B - D   2.154936 0.7177         -11.111791 15.4216641
B - E   9.061994 0.1539          -4.204734 22.3287219
C - D  14.438186 0.0364       *   1.171459 27.7049140
C - E  21.345244 0.0060      **   8.078516 34.6119718
D - E   6.907058 0.2643          -6.359670 20.1737855

Diseño completamente al azar en arreglo factorial con tres vías

Planeamiento

Crear un libro de campo con el paquete agricolae

trt <- c(2,2,2)

r <- 5

salida <- agricolae::design.ab(trt = trt,
                               r = r,
                               design = "crd",
                               serie = 3,
                               seed = 123,
                               kinds = "Super-Duper",
                               randomization = TRUE)
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Guardar el libro generado

write.table(salida$book,
            "books/crd3w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida$book,
           "books/crd3w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
agricolaeplotr::plot_design.factorial_crd(
  design = salida,
  factor_name = "A",
  ncols = 8,
  nrows = 5,
  reverse_y = TRUE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas")
agricolaeplotr::plot_design.factorial_crd(
  design = salida,
  factor_name = "B",
  ncols = 8,
  nrows = 5,
  reverse_y = TRUE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Columnas",
       y = "Filas")
agricolaeplotr::plot_design.factorial_crd(
  design = salida,
  factor_name = "C",
  ncols = 8,
  nrows = 5,
  reverse_y = TRUE) +
  labs(fill = "Dosis de Potasio",
       x = "Columnas",
       y = "Filas")

Crear un libro de campo con el paquete edibble

menu_factorial()
design("Factorial Design") %>%
  set_units(unit = 16) %>%
  set_trts(trt1 = 4,
           trt2 = 2) %>%
  allot_trts(~unit) %>%
  assign_trts("random", seed = 903) %>%
  serve_table()
crd <- takeout(menu_factorial(trt = c(2,2,2),
                              r = 5,
                              seed = 123))
crd %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
crd2 <- design("Factorial Design") %>%
  set_units(unit = 40) %>%
  set_trts(trt1 = 2,
           trt2 = 2,
           trt3 = 2) %>%
  allot_trts(~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
crd2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
deggust::autoplot(crd)
plot(crd)

Análisis de DCA en arreglo factorial con tres vías

Importación de datos

archivos <- list.files(pattern = "datos dca 3w.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor) %>%
  rename("rep" = "bloque")

Creación del modelo lineal

modelo.dca1 <- lm(respuesta ~ a*b*c, data = data)
modelo.dca2 <- lm(respuesta ~ a*b+c, data = data)
modelo.dca3 <- lm(respuesta ~ a*c+b, data = data)
modelo.dca4 <- lm(respuesta ~ a+b*c, data = data)
modelo.dca5 <- lm(respuesta ~ a+b+c, data = data)
broom::glance(modelo.dca1) %>%
  bind_rows(broom::glance(modelo.dca2),
            broom::glance(modelo.dca3),
            broom::glance(modelo.dca4),
            broom::glance(modelo.dca5)) %>%
  dplyr::mutate(Modelo = c("a*b*c", "a*b+c", "a*c+b",
                    "a+b*c", "a+b+c")) %>%
  dplyr::relocate(Modelo) %>%
  dplyr::select(Modelo, AIC, BIC) %>%
  dplyr::arrange(BIC) %>%
  gt()
Modelo AIC BIC
a+b+c 143.9536 149.8439
a+b*c 143.4107 150.4790
a*c+b 145.0689 152.1372
a*b+c 145.5644 152.6328
a*b*c 146.0778 156.6803
 
modelo.dca <- lm(respuesta ~ a+b+c, data = data)

Definición del modelo

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2\]

summary(modelo.dca)

Call:
lm(formula = respuesta ~ a + b + c, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-6.500 -3.125 -0.250  3.625  7.500 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   14.000      1.763   7.941 1.31e-07 ***
aa2            3.500      1.763   1.985  0.06101 .  
bb2            6.500      1.763   3.687  0.00146 ** 
cc2            2.000      1.763   1.134  0.27004    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.319 on 20 degrees of freedom
Multiple R-squared:  0.4848,    Adjusted R-squared:  0.4075 
F-statistic: 6.273 on 3 and 20 DF,  p-value: 0.003551

Verificación visual de los supuestos del modelo

performance::check_model(modelo.dca)
ggResidpanel::resid_panel(modelo.dca)
influence.measures(modelo.dca)
Influence measures of
     lm(formula = respuesta ~ a + b + c, data = data) :

      dfb.1_ dfb.aa2 dfb.bb2 dfb.cc2     dffit cov.r   cook.d   hat inf
1  -2.27e-01  0.2271  0.2271  0.2271  4.54e-01 1.193 5.15e-02 0.167    
2  -8.32e-02 -0.0832  0.0832  0.0832 -1.66e-01 1.431 7.24e-03 0.167    
3   1.97e-01 -0.1974  0.1974 -0.1974  3.95e-01 1.254 3.94e-02 0.167    
4   1.68e-01 -0.1683 -0.1683  0.1683  3.37e-01 1.310 2.90e-02 0.167    
5   7.72e-17  0.2271  0.2271 -0.2271  4.54e-01 1.193 5.15e-02 0.167    
6  -8.05e-17 -0.1974  0.1974 -0.1974 -3.95e-01 1.254 3.94e-02 0.167    
7   1.22e-16 -0.2573  0.2573  0.2573  5.15e-01 1.125 6.51e-02 0.167    
8  -7.06e-01  0.3528  0.3528  0.3528 -7.06e-01 0.900 1.16e-01 0.167    
9   1.68e-01 -0.1683 -0.1683 -0.1683 -3.37e-01 1.310 2.90e-02 0.167    
10 -8.32e-02 -0.0832  0.0832  0.0832 -1.66e-01 1.431 7.24e-03 0.167    
11  4.58e-01 -0.4582  0.4582 -0.4582  9.16e-01 0.663 1.81e-01 0.167    
12  0.00e+00  0.0000  0.0000  0.0000  8.59e-17 1.473 1.94e-33 0.167    
13 -1.20e-16 -0.3528 -0.3528  0.3528 -7.06e-01 0.900 1.16e-01 0.167    
14  1.31e-16  0.3200 -0.3200  0.3200  6.40e-01 0.978 9.73e-02 0.167    
15 -1.84e-16  0.3866 -0.3866 -0.3866 -7.73e-01 0.821 1.36e-01 0.167    
16 -4.54e-01  0.2271  0.2271  0.2271 -4.54e-01 1.193 5.15e-02 0.167    
17  5.54e-02 -0.0554 -0.0554 -0.0554 -1.11e-01 1.454 3.22e-03 0.167    
18  8.32e-02  0.0832 -0.0832 -0.0832  1.66e-01 1.431 7.24e-03 0.167    
19 -2.77e-02  0.0277 -0.0277  0.0277 -5.53e-02 1.469 8.04e-04 0.167    
20  0.00e+00  0.0000  0.0000  0.0000  8.59e-17 1.473 1.94e-33 0.167    
21 -1.88e-17 -0.0554 -0.0554  0.0554 -1.11e-01 1.454 3.22e-03 0.167    
22  5.70e-17  0.1396 -0.1396  0.1396  2.79e-01 1.358 2.01e-02 0.167    
23 -1.52e-16  0.3200 -0.3200 -0.3200 -6.40e-01 0.978 9.73e-02 0.167    
24  4.54e-01 -0.2271 -0.2271 -0.2271  4.54e-01 1.193 5.15e-02 0.167    
influenceIndexPlot(modelo.dca)

Cumplimiento de supuestos del modelo lineal general

Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dca,
                 reps = 5000,
                 max.lag = 5)
 lag Autocorrelation D-W Statistic p-value
   1     -0.30428954      2.522788  0.1232
   2      0.03753351      1.752011  0.6948
   3     -0.01005362      1.797587  0.9592
   4      0.08579088      1.579088  0.5168
   5     -0.27949062      2.266756  0.0960
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dca, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dca
DW = 2.5228, p-value = 0.1185
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dca))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dca)
W = 0.96315, p-value = 0.5049
ad.test(rstudent(modelo.dca))

    Anderson-Darling normality test

data:  rstudent(modelo.dca)
A = 0.31136, p-value = 0.5281
lillie.test(rstudent(modelo.dca))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dca)
D = 0.099285, p-value = 0.7825
ks.test(rstudent(modelo.dca), "pnorm",
        alternative = "two.sided")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dca)
D = 0.10882, p-value = 0.9388
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dca))

    Cramer-von Mises normality test

data:  rstudent(modelo.dca)
W = 0.044742, p-value = 0.5783
pearson.test(rstudent(modelo.dca))

    Pearson chi-square normality test

data:  rstudent(modelo.dca)
P = 4.6667, p-value = 0.4579
sf.test(rstudent(modelo.dca))

    Shapiro-Francia normality test

data:  rstudent(modelo.dca)
W = 0.97286, p-value = 0.6459

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento

ncvTest(modelo.dca)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.1309303, Df = 1, p = 0.71747
bptest(modelo.dca)

    studentized Breusch-Pagan test

data:  modelo.dca
BP = 4.1399, df = 3, p-value = 0.2467
bptest(modelo.dca, studentize = F)

    Breusch-Pagan test

data:  modelo.dca
BP = 1.9734, df = 3, p-value = 0.5779
olsrr::ols_test_breusch_pagan(modelo.dca)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

                Data                  
 -------------------------------------
 Response : respuesta 
 Variables: fitted values of respuesta 

        Test Summary         
 ----------------------------
 DF            =    1 
 Chi2          =    0.1309303 
 Prob > Chi2   =    0.7174694

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que no se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dca %>% gvlma()

Call:
lm(formula = respuesta ~ a + b + c, data = data)

Coefficients:
(Intercept)          aa2          bb2          cc2  
       14.0          3.5          6.5          2.0  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                      Value p-value                Decision
Global Stat        2.211805  0.6969 Assumptions acceptable.
Skewness           0.005698  0.9398 Assumptions acceptable.
Kurtosis           1.095441  0.2953 Assumptions acceptable.
Link Function      1.051042  0.3053 Assumptions acceptable.
Heteroscedasticity 0.059624  0.8071 Assumptions acceptable.

Análisis de varianza

\[Y_{ijk} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \epsilon_{ijk}\]

\[\hat{Y}_{ijk} = \mu + \tau_{i} + \beta_{j} + \gamma_{k}\]

Para el factor A (Variedad):

\(H_0: \tau_{A1} = \tau_{A2} = 0\)

\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)

Para el factor B (Densidad):

\(H_0: \beta_{B1} = \beta_{B2} = 0\)

\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_j \neq 0\text{; en al menos un nivel del factor B.}\)

Para el factor C (Época):

\(H_0: \gamma_{C1} = \gamma_{C2} = 0\)

\(H_1: \text{En al menos un nivel del factor C el } \gamma \text{ es diferente a los demás.}\)

\(H_1: \gamma_k \neq 0\text{; en al menos un nivel del factor C.}\)

anova(modelo.dca, test = "F")
Analysis of Variance Table

Response: respuesta
          Df Sum Sq Mean Sq F value   Pr(>F)   
a          1   73.5   73.50  3.9410 0.061007 . 
b          1  253.5  253.50 13.5925 0.001461 **
c          1   24.0   24.00  1.2869 0.270041   
Residuals 20  373.0   18.65                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 1, 20)
[1] 4.351244

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 1, 20)
[1] 4.351244

Valor de la tabla de F para el factor C con una significancia de 0.05.

qf(0.95, 1, 20)
[1] 4.351244
broom::augment(modelo.dca)
# A tibble: 24 × 10
   respuesta a     b     c     .fitted .resid  .hat .sigma .cooksd .std.resid
       <dbl> <fct> <fct> <fct>   <dbl>  <dbl> <dbl>  <dbl>   <dbl>      <dbl>
 1        30 a2    b2    c2       26      4   0.167   4.32 0.0515       1.01 
 2        16 a2    b1    c1       17.5   -1.5 0.167   4.41 0.00724     -0.380
 3        24 a1    b2    c1       20.5    3.5 0.167   4.34 0.0394       0.888
 4        19 a1    b1    c2       16      3   0.167   4.37 0.0290       0.761
 5        28 a2    b2    c1       24      4   0.167   4.32 0.0515       1.01 
 6        16 a2    b1    c2       19.5   -3.5 0.167   4.34 0.0394      -0.888
 7        27 a1    b2    c2       22.5    4.5 0.167   4.28 0.0651       1.14 
 8         8 a1    b1    c1       14     -6   0.167   4.17 0.116       -1.52 
 9        23 a2    b2    c2       26     -3   0.167   4.37 0.0290      -0.761
10        16 a2    b1    c1       17.5   -1.5 0.167   4.41 0.00724     -0.380
# ℹ 14 more rows
broom::tidy(modelo.dca)
# A tibble: 4 × 5
  term        estimate std.error statistic     p.value
  <chr>          <dbl>     <dbl>     <dbl>       <dbl>
1 (Intercept)     14        1.76      7.94 0.000000131
2 aa2              3.5      1.76      1.99 0.0610     
3 bb2              6.5      1.76      3.69 0.00146    
4 cc2              2        1.76      1.13 0.270      
anova(modelo.dca, test = "F") %>% 
  as.data.frame() %>%
  dplyr::transmute(porc.varexp = `Sum Sq`/sum(`Sum Sq`)*100)
          porc.varexp
a           10.151934
b           35.013812
c            3.314917
Residuals   51.519337
shadeDist(qf(0.95, 1, 20), "df",1, 20, lower.tail = F)

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los efectos de los niveles del factor A son estadísticamente similares a 0.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor B.

Con respecto al Factor C: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los efectos de los niveles del factor C son estadísticamente similares a 0.

agricolae::cv.model(modelo.dca)
[1] 21.59282

Comparaciones de medias

Para los niveles del factor B:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

Prueba de LSD

Para los niveles del Factor B:

agricolae::LSD.test(modelo.dca, trt = "b", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dca ~ "b"

LSD t Test for respuesta 

Mean Square Error:  18.65 

b,  means and individual ( 95 %) CI

   respuesta      std  r      LCL      UCL Min Max
b1     16.75 4.594958 12 14.14951 19.35049   8  25
b2     23.25 4.653933 12 20.64951 25.85049  16  30

Alpha: 0.05 ; DF Error: 20
Critical Value of t: 2.085963 

least Significant Difference: 3.677651 

Treatments with the same letter are not significantly different.

   respuesta groups
b2     23.25      a
b1     16.75      b
agricolae::LSD.test(modelo.dca, trt = "b", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dca ~ "b"

LSD t Test for respuesta 

Mean Square Error:  18.65 

b,  means and individual ( 95 %) CI

   respuesta      std  r      LCL      UCL Min Max
b1     16.75 4.594958 12 14.14951 19.35049   8  25
b2     23.25 4.653933 12 20.64951 25.85049  16  30

Alpha: 0.05 ; DF Error: 20
Critical Value of t: 2.085963 

Comparison between treatments means

        difference pvalue signif.       LCL       UCL
b1 - b2       -6.5 0.0015      ** -10.17765 -2.822349

DCA factorial pxqxr con el paquete ExpDes

attach(data)
fat3.crd(factor1 = a,
         factor2 =  b,
         factor3 =  c,
         resp =  respuesta,
         quali = c(TRUE,
                   TRUE,
                   TRUE),
         mcomp = "tukey", 
         fac.names = c("Variedades",
                       "Densidad",
                       "Época"),
         sigT = 0.05,
         sigF = 0.05,
         unfold = NULL)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Variedades 
FACTOR 2:  Densidad 
FACTOR 3:  Época 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                          DF    SS    MS      Fc  Pr>Fc
Variedades                 1  73.5  73.5  4.0274  0.062
Densidad                   1 253.5 253.5 13.8904 0.0018
Época                      1  24.0    24  1.3151 0.2683
Variedades*Densidad        1   6.0     6  0.3288 0.5744
Variedades*Época           1  13.5  13.5  0.7397 0.4025
Densidad*Época             1  37.5  37.5  2.0548  0.171
Variedades*Densidad*Época  1  24.0    24  1.3151 0.2683
Residuals                 16 292.0 18.25               
Total                     23 724.0                     
------------------------------------------------------------------------
CV = 21.36 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.1816204 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------

No significant interaction: analyzing the simple effect
------------------------------------------------------------------------
Variedades
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels Means
1     a1 18.25
2     a2 21.75
------------------------------------------------------------------------
Densidad
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    b2      23.25 
 b   b1      16.75 
------------------------------------------------------------------------

Época
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels Means
1     c1    19
2     c2    21
------------------------------------------------------------------------

DCA factorial pxqxr con tratamientos adicionales con el paquete ExpDes

data(ex6)
attach(ex6)
data(respAd)
fat3.ad.crd(
  factor1 = fatorA,
  factor2 =  fatorB,
  factor3 =  fatorC,
  repet =  rep,
  resp =  resp,
  respAd =  respAd,
  quali = c(TRUE,TRUE, TRUE),
  mcomp = "tukey",
  fac.names =
    c("Factor A",
      "Factor B",
      "Factor C"), 
  sigT = 0.05, 
  sigF = 0.05)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Factor A 
FACTOR 2:  Factor B 
FACTOR 3:  Factor C 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                           DF       SS      MS     Fc  Pr>Fc
Factor A                    1  0.51042 0.51042 0.9381 0.3456
Factor B                    1  0.51042 0.51042 0.9381 0.3456
Factor C                    1  0.05042 0.05042 0.0927 0.7643
Factor A*Factor B           1  2.87042 2.87042 5.2758 0.0338
Factor A*Factor C           1  0.12042 0.12042 0.2213 0.6437
Factor B*Factor C           1  0.40042 0.40042  0.736 0.4022
Factor A*Factor B*Factor C  1  0.02042 0.02042 0.0375 0.8486
Ad vs Factorial             1  0.57042 0.57042 1.0484 0.3194
Residuals                  18  9.79333 0.54407              
Total                      26 14.84667                      
------------------------------------------------------------------------
CV = 7.29 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.768125 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
According to the F test, the means of the two groups are statistical equal.
              Means
Additional 10.53333
Factorial  10.07083
------------------------------------------------------------------------



Significant Factor A*Factor B  interaction: analyzing the interaction
------------------------------------------------------------------------

Analyzing  Factor A  inside of each level of  Factor B 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                    DF      SS      MS     Fc Pr>Fc
Factor A:Factor B 1  1 2.90083 2.90083 5.3317 0.033
Factor A:Factor B 2  1 0.48000 0.48000 0.8822  0.36
Residuals           18 9.79333 0.54407             
------------------------------------------------------------------------



 Factor A  inside of the level  1  of  Factor B 
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    1   10.41667 
 b   2   9.433333 
------------------------------------------------------------------------


 Factor A  inside of the level  2  of  Factor B 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels    Means
1      1 10.01667
2      2 10.41667
------------------------------------------------------------------------



Analyzing  Factor B  inside of each level of  Factor A 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                    DF      SS      MS     Fc Pr>Fc
Factor B:Factor A 1  1 0.48000 0.48000 0.8822  0.36
Factor B:Factor A 2  1 2.90083 2.90083 5.3317 0.033
Residuals           18 9.79333 0.54407             
------------------------------------------------------------------------



 Factor B  inside of the level  1  of  Factor A 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels    Means
1      1 10.41667
2      2 10.01667
------------------------------------------------------------------------


 Factor B  inside of the level  2  of  Factor A 
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    2   10.41667 
 b   1   9.433333 
------------------------------------------------------------------------

Analizing the effect of the factor  Factor C 
------------------------------------------------------------------------
Factor C
According to the F test, the means of this factor are not different.
------------------------------------------------------------------------
  Niveis   Medias
1      1 10.02500
2      2 10.11667
------------------------------------------------------------------------

Prueba de Dunnet

data <- ex6 %>%
  mutate_if(is.integer, factor)%>% 
  dplyr::select(fatorA, fatorB, fatorC, rep, resp) %>%
  dplyr::bind_rows(data.frame(fatorA = "ad",
                       fatorB = "ad",
                       fatorC = "ad",
                       rep = factor(c(1:3)),
                       resp = respAd)) %>%
  dplyr::mutate(tratamiento =
                  interaction(fatorA,
                              fatorB,
                              fatorC))
modelo <- lm(resp ~ tratamiento, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = resp ~ tratamiento, data = data)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)
2.1.1 - 1.1.1 == 0    -0.90000    0.60226  -1.494    0.564
1.2.1 - 1.1.1 == 0    -0.20000    0.60226  -0.332    1.000
2.2.1 - 1.1.1 == 0     0.40000    0.60226   0.664    0.984
1.1.2 - 1.1.1 == 0     0.43333    0.60226   0.720    0.975
2.1.2 - 1.1.1 == 0    -0.63333    0.60226  -1.052    0.851
1.2.2 - 1.1.1 == 0    -0.16667    0.60226  -0.277    1.000
2.2.2 - 1.1.1 == 0     0.03333    0.60226   0.055    1.000
ad.ad.ad - 1.1.1 == 0  0.33333    0.60226   0.553    0.995
(Adjusted p values reported -- single-step method)

Especificar el nivel del tratamiento adicional

data$tratamiento <-
  relevel(data$tratamiento,"ad.ad.ad")
modelo <- lm(resp ~ tratamiento, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = resp ~ tratamiento, data = data)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)
1.1.1 - ad.ad.ad == 0 -0.33333    0.60226  -0.553    0.995
2.1.1 - ad.ad.ad == 0 -1.23333    0.60226  -2.048    0.256
1.2.1 - ad.ad.ad == 0 -0.53333    0.60226  -0.886    0.928
2.2.1 - ad.ad.ad == 0  0.06667    0.60226   0.111    1.000
1.1.2 - ad.ad.ad == 0  0.10000    0.60226   0.166    1.000
2.1.2 - ad.ad.ad == 0 -0.96667    0.60226  -1.605    0.491
1.2.2 - ad.ad.ad == 0 -0.50000    0.60226  -0.830    0.947
2.2.2 - ad.ad.ad == 0 -0.30000    0.60226  -0.498    0.997
(Adjusted p values reported -- single-step method)

Diseño de bloques completos al azar en arreglo factorial con dos vías

Planeamiento

Crear un libro de campo con el paquete agricolae

trt <- c(2,2,2)

r <- 5

salida <- agricolae::design.ab(trt = trt,
                               r = r,
                               design = "rcbd",
                               serie = 3,
                               seed = 123,
                               kinds = "Mersenne-Twister",
                               randomization = TRUE)
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Guardar el libro generado

write.table(salida$book,
            "books/rcbd3w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida$book,
           "books/rcbd3w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)

Guardar el libro de campo generado

write.table(salida %>% zigzag(),
            "books/rcbd3w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida %>% zigzag(),
           "books/rcbd3w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
# Función no aplicable debido a bugs

agricolaeplotr::plot_design.factorial_rcbd(
  design = salida,
  factor_name = "A",
  reverse_y = TRUE,
  reverse_x = TRUE) +
  labs(fill = "Variedades",
       x = "Bloque",
       y = "Columna")
plots <- as.numeric(salida$book[, 1])
mat <- matrix(plots, byrow = TRUE, ncol = length(levels(interaction(salida$book$A,salida$book$B,salida$book$C))))
dims <- dim(mat)
table <- as.table(mat)
rownames(table) <- 1:dims[1]
colnames(table) <- 1:dims[2]
table <- as.data.frame(salida$book) %>%
  rename("Bloque" = "block") %>%
  dplyr::mutate(Columna = rep(c(1:length(levels(interaction(salida$book$A,salida$book$B,salida$book$C)))),5))

ggplot(table, aes(x = Columna, y = Bloque)) + 
      geom_tile(aes(fill = A)) + 
      theme_bw() + theme(line = element_blank()) + geom_text(aes_string(label = "plots"), 
      colour = "black")
plot_design.factorial_rcbd_corregido(design = salida,
  factor_name = "A",
  reverse_y = TRUE,
  reverse_x = TRUE) +
  labs(fill = "Variedades",
       x = "Columna",
       y = "Bloque")
# Función no aplicable debido a bugs

agricolaeplotr::plot_design.factorial_rcbd(
  design = salida,
  factor_name = "B",
  reverse_y = FALSE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Bloque",
       y = "Columnas")
plot_design.factorial_rcbd_corregido(design = salida,
  factor_name = "B",
  reverse_y = TRUE,
  reverse_x = TRUE) +
  labs(fill = "Densidad",
       x = "Columna",
       y = "Bloque")
plot_design.factorial_rcbd_corregido(design = salida,
  factor_name = "C",
  reverse_y = TRUE,
  reverse_x = TRUE) +
  labs(fill = "Época",
       x = "Columna",
       y = "Bloque")

Crear un libro de campo con el paquete edibble

menu_factorial()
design("Factorial Design") %>%
  set_units(unit = 162) %>%
  set_trts(trt1 = 6,
           trt2 = 9) %>%
  allot_trts(~unit) %>%
  assign_trts("random", seed = 865) %>%
  serve_table()
rcbd <- takeout(menu_factorial(trt = c(2,2,2),
                              r = 5 ,
                              seed = 123))
rcbd %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
rcbd2 <- design("Factorial Design") %>%
  set_units(block = 5,
            unit = nested_in(block, 8)) %>%
  set_trts(trt1 = 2,
           trt2 = 2,
           trt3 = 2) %>%
  allot_trts(trt1 + trt2 + trt3 ~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
rcbd2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
write.table(rcbd2 %>% as.data.frame(),
            "books/rcbd3w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(rcbd2 %>% as.data.frame(),
           "books/rcbd3w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
deggust::autoplot(rcbd2)
plot(rcbd2)

Análisis de DBCA en arreglo factorial con tres vías

Importación de datos

archivos <- list.files(pattern = "datos dbca 3w.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor) %>%
  mutate(bloque = factor(bloque))

Creación del modelo lineal

modelo.dbca1 <- lm(respuesta ~ a*b*c+bloque, data = data)
modelo.dbca2 <- lm(respuesta ~ a*b+c+bloque, data = data)
modelo.dbca3 <- lm(respuesta ~ a*c+b+bloque, data = data)
modelo.dbca4 <- lm(respuesta ~ a+b*c+bloque, data = data)
modelo.dbca5 <- lm(respuesta ~ a+b+c+bloque, data = data)
broom::glance(modelo.dbca1) %>%
  bind_rows(broom::glance(modelo.dbca2),
            broom::glance(modelo.dbca3),
            broom::glance(modelo.dbca4),
            broom::glance(modelo.dbca5)) %>%
  dplyr::mutate(Modelo = c("a*b*c", "a*b+c", "a*c+b",
                    "a+b*c", "a+b+c")) %>%
  dplyr::select(Modelo, AIC, BIC) %>%
  dplyr::arrange(BIC) %>%
  dplyr::mutate(Mérito = 1:n()) %>%
  dplyr::relocate(Mérito, Modelo) %>%
  gt()
Mérito Modelo AIC BIC
1 a+b+c 146.9014 155.1478
2 a+b*c 146.2379 155.6624
3 a*c+b 147.9762 157.4007
4 a*b+c 148.4946 157.9191
5 a*b*c 148.7254 161.6840
 
modelo.dbca <- lm(respuesta ~ a + b + c + bloque, data = data)

Definición del modelo

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2 + \beta_4*Bloq_{II} + \beta_5*Bloq_{III} + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2 + \beta_4*Bloq_{II} + \beta_5*Bloq_{III}\]

summary(modelo.dbca)

Call:
lm(formula = respuesta ~ a + b + c + bloque, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-7.000 -2.625 -0.250  2.625  8.500 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   15.000      2.227   6.736 2.59e-06 ***
aa2            3.500      1.818   1.925  0.07017 .  
bb2            6.500      1.818   3.575  0.00216 ** 
cc2            2.000      1.818   1.100  0.28581    
bloque2       -2.000      2.227  -0.898  0.38095    
bloque3       -1.000      2.227  -0.449  0.65873    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.453 on 18 degrees of freedom
Multiple R-squared:  0.5069,    Adjusted R-squared:  0.3699 
F-statistic: 3.701 on 5 and 18 DF,  p-value: 0.01768

Verificación visual de los supuestos del modelo

performance::check_model(modelo.dbca)
ggResidpanel::resid_panel(modelo.dbca)
influence.measures(modelo.dbca)
Influence measures of
     lm(formula = respuesta ~ a + b + c + bloque, data = data) :

      dfb.1_ dfb.aa2 dfb.bb2 dfb.cc2  dfb.blq2  dfb.blq3     dffit cov.r
1  -4.93e-16  0.1812  0.1812  0.1812 -2.22e-01 -2.22e-01  4.44e-01 1.530
2  -2.45e-01 -0.1502  0.1502  0.1502  1.84e-01  1.84e-01 -3.68e-01 1.631
3   2.45e-01 -0.1502  0.1502 -0.1502 -1.84e-01 -1.84e-01  3.68e-01 1.631
4   1.95e-01 -0.1197 -0.1197  0.1197 -1.47e-01 -1.47e-01  2.93e-01 1.717
5   1.48e-01  0.1812  0.1812 -0.1812 -2.22e-01 -2.22e-01  4.44e-01 1.530
6  -2.27e-01 -0.2780  0.2780 -0.2780  3.40e-01  3.40e-01 -6.81e-01 1.172
7   1.74e-01 -0.2128  0.2128  0.2128 -2.61e-01 -2.61e-01  5.21e-01 1.419
8  -1.13e+00  0.4600  0.4600  0.4600  5.63e-01  5.63e-01 -1.13e+00 0.559
9   1.47e-01 -0.1197 -0.1197 -0.1197 -1.47e-01 -4.48e-17 -2.93e-01 1.717
10 -1.21e-02 -0.0297  0.0297  0.0297 -3.64e-02  1.01e-18 -7.28e-02 1.868
11  2.41e-01 -0.5908  0.5908 -0.5908  7.24e-01  2.01e-17  1.45e+00 0.285
12  2.43e-02 -0.0595 -0.0595  0.0595  7.29e-02  1.21e-17  1.46e-01 1.837
13  1.27e-01 -0.3119 -0.3119  0.3119 -3.82e-01  5.30e-17 -7.64e-01 1.043
14 -1.72e-01  0.4207 -0.4207  0.4207  5.15e-01  1.00e-16  1.03e+00 0.670
15  1.42e-01  0.3468 -0.3468 -0.3468 -4.25e-01 -1.06e-16 -8.50e-01 0.915
16 -2.22e-01  0.1812  0.1812  0.1812 -2.22e-01  1.23e-17 -4.44e-01 1.530
17  7.29e-02 -0.0595 -0.0595 -0.0595 -1.40e-17 -7.29e-02 -1.46e-01 1.837
18  3.65e-02  0.0895 -0.0895 -0.0895 -5.27e-17  1.10e-01  2.19e-01 1.786
19 -1.21e-02  0.0297 -0.0297  0.0297  1.05e-17 -3.64e-02 -7.28e-02 1.868
20  0.00e+00  0.0000  0.0000  0.0000  0.00e+00  0.00e+00  9.69e-17 1.879
21  2.43e-02 -0.0595 -0.0595  0.0595  1.40e-17 -7.29e-02 -1.46e-01 1.837
22 -6.13e-02  0.1502 -0.1502  0.1502 -3.54e-17  1.84e-01  3.68e-01 1.631
23  1.42e-01  0.3468 -0.3468 -0.3468  0.00e+00 -4.25e-01 -8.50e-01 0.915
24  3.00e-01 -0.2450 -0.2450 -0.2450 -8.66e-17  3.00e-01  6.00e-01 1.298
     cook.d  hat inf
1  3.36e-02 0.25    
2  2.33e-02 0.25    
3  2.33e-02 0.25    
4  1.49e-02 0.25    
5  3.36e-02 0.25    
6  7.56e-02 0.25    
7  4.58e-02 0.25    
8  1.83e-01 0.25   *
9  1.49e-02 0.25    
10 9.34e-04 0.25    
11 2.70e-01 0.25    
12 3.73e-03 0.25    
13 9.34e-02 0.25    
14 1.58e-01 0.25    
15 1.13e-01 0.25    
16 3.36e-02 0.25    
17 3.73e-03 0.25    
18 8.40e-03 0.25    
19 9.34e-04 0.25    
20 1.66e-33 0.25    
21 3.73e-03 0.25    
22 2.33e-02 0.25    
23 1.13e-01 0.25    
24 5.98e-02 0.25    
influenceIndexPlot(modelo.dbca)

Cumplimiento de supuestos del modelo lineal general

Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dbca,
                 reps = 5000,
                 max.lag = 5)
 lag Autocorrelation D-W Statistic p-value
   1     -0.36834734      2.666667  0.1400
   2      0.03221289      1.763305  0.4492
   3     -0.09453782      1.981793  0.8552
   4      0.04481793      1.689076  0.5556
   5     -0.28641457      2.326331  0.1056
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dbca, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dbca
DW = 2.6667, p-value = 0.1332
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dbca))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dbca)
W = 0.9787, p-value = 0.8708
ad.test(rstudent(modelo.dbca))

    Anderson-Darling normality test

data:  rstudent(modelo.dbca)
A = 0.20956, p-value = 0.8435
lillie.test(rstudent(modelo.dbca))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dbca)
D = 0.084326, p-value = 0.93
ks.test(rstudent(modelo.dbca), "pnorm",
        alternative = "two.sided")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dbca)
D = 0.089208, p-value = 0.991
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))

    Cramer-von Mises normality test

data:  rstudent(modelo.dbca)
W = 0.027358, p-value = 0.8739
pearson.test(rstudent(modelo.dbca))

    Pearson chi-square normality test

data:  rstudent(modelo.dbca)
P = 4.6667, p-value = 0.4579
sf.test(rstudent(modelo.dbca))

    Shapiro-Francia normality test

data:  rstudent(modelo.dbca)
W = 0.97797, p-value = 0.7739

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento

ncvTest(modelo.dbca)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.1259105, Df = 1, p = 0.72271
bptest(modelo.dbca)

    studentized Breusch-Pagan test

data:  modelo.dbca
BP = 4.9286, df = 5, p-value = 0.4247
bptest(modelo.dbca, studentize = F)

    Breusch-Pagan test

data:  modelo.dbca
BP = 3.6045, df = 5, p-value = 0.6076
olsrr::ols_test_breusch_pagan(modelo.dbca)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

                Data                  
 -------------------------------------
 Response : respuesta 
 Variables: fitted values of respuesta 

        Test Summary         
 ----------------------------
 DF            =    1 
 Chi2          =    0.1259105 
 Prob > Chi2   =    0.7227104

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que no se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dbca %>% gvlma()

Call:
lm(formula = respuesta ~ a + b + c + bloque, data = data)

Coefficients:
(Intercept)          aa2          bb2          cc2      bloque2      bloque3  
       15.0          3.5          6.5          2.0         -2.0         -1.0  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                     Value p-value                Decision
Global Stat        0.38057  0.9840 Assumptions acceptable.
Skewness           0.04653  0.8292 Assumptions acceptable.
Kurtosis           0.28871  0.5910 Assumptions acceptable.
Link Function      0.02489  0.8747 Assumptions acceptable.
Heteroscedasticity 0.02044  0.8863 Assumptions acceptable.

Análisis de varianza

\[Y_{ijkl} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_l + \epsilon_{ijkl}\]

\[\hat{Y}_{ijkl} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_l\]

Dónde:

\(Y_{ijkl}\) = Valor observado de la variable respuesta.

\(\hat{Y}_{ijkl}\) = Valor ajustado de la variable respuesta.

\(\mu\) = Promedio observado de la variable respuesta.

\(\tau_{i}\) = Efecto del i-ésimo nivel del factor A.

\(\beta_{j}\) = Efecto del j-ésimo nivel del factor B.

\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor C.

\(\delta_{l}\) = Efecto del l-ésimo nivel del factor Bloque.

\(\epsilon_{ijkl}\) = Residuo observado del modelo.

Pruebas de hipótesis

Para el factor A (Variedad):

\(H_0: \tau_{A1} = \tau_{A2} = 0\)

\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_{i} \neq 0\text{; en al menos un nivel del factor A.}\)

Para el factor B (Densidad):

\(H_0: \beta_{B1} = \beta_{B2} = 0\)

\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_{j} \neq 0\text{; en al menos un nivel del factor B.}\)

Para el factor C (Época):

\(H_0: \gamma_{C1} = \gamma_{C2} = 0\)

\(H_1: \text{En al menos un nivel del factor C el } \gamma \text{ es diferente a los demás.}\)

\(H_1: \gamma_{k} \neq 0\text{; en al menos un nivel del factor C.}\)

anova(modelo.dbca, test = "F")
Analysis of Variance Table

Response: respuesta
          Df Sum Sq Mean Sq F value   Pr(>F)   
a          1   73.5  73.500  3.7059 0.070172 . 
b          1  253.5 253.500 12.7815 0.002164 **
c          1   24.0  24.000  1.2101 0.285810   
bloque     2   16.0   8.000  0.4034 0.673961   
Residuals 18  357.0  19.833                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 1, 18)
[1] 4.413873

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 1, 18)
[1] 4.413873

Valor de la tabla de F para el factor C con una significancia de 0.05.

qf(0.95, 1, 18)
[1] 4.413873

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor A tienen un efecto sobre el rendimiento estadísticamente similar a 0.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor B.

Con respecto al Factor C: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor C tienen un efecto sobre el rendimiento estadísticamente similar a 0.

agricolae::cv.model(modelo.dbca)
[1] 22.26732

Comparaciones de medias

Para los niveles del factor B:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

Prueba de LSD

Para los niveles del Factor B:

agricolae::LSD.test(modelo.dbca, trt = "b", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dbca ~ "b"

LSD t Test for respuesta 

Mean Square Error:  19.83333 

b,  means and individual ( 95 %) CI

   respuesta      std  r      LCL      UCL Min Max
b1     16.75 4.594958 12 14.04905 19.45095   8  25
b2     23.25 4.653933 12 20.54905 25.95095  16  30

Alpha: 0.05 ; DF Error: 18
Critical Value of t: 2.100922 

least Significant Difference: 3.819726 

Treatments with the same letter are not significantly different.

   respuesta groups
b2     23.25      a
b1     16.75      b
agricolae::LSD.test(modelo.dbca, trt = "b", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dbca ~ "b"

LSD t Test for respuesta 

Mean Square Error:  19.83333 

b,  means and individual ( 95 %) CI

   respuesta      std  r      LCL      UCL Min Max
b1     16.75 4.594958 12 14.04905 19.45095   8  25
b2     23.25 4.653933 12 20.54905 25.95095  16  30

Alpha: 0.05 ; DF Error: 18
Critical Value of t: 2.100922 

Comparison between treatments means

        difference pvalue signif.       LCL       UCL
b1 - b2       -6.5 0.0022      ** -10.31973 -2.680274

Análisis de varianza de un DBCA factorial con tres vías con el paquete ExpDes

attach(data)
fat3.rbd(factor1 = a,
         factor2 = b,
         factor3 = c,
         block = bloque,
         resp =  respuesta,
         quali = c(TRUE,
                   TRUE,
                   TRUE),
         mcomp = "lsd", 
         fac.names = c("Variedades",
                       "Densidad",
                       "Época"),
         sigT = 0.05,
         sigF = 0.05)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Variedades 
FACTOR 2:  Densidad 
FACTOR 3:  Época 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                          DF    SS       MS      Fc  Pr>Fc
Block                      2  16.0        8  0.4058  0.674
Variedades                 1  73.5     73.5  3.7283  0.074
Densidad                   1 253.5    253.5 12.8587  0.003
Época                      1  24.0       24  1.2174 0.2885
Variedades*Densidad        1   6.0        6  0.3043 0.5899
Variedades*Época           1  13.5     13.5  0.6848 0.4218
Densidad*Época             1  37.5     37.5  1.9022 0.1895
Variedades*Densidad*Época  1  24.0       24  1.2174 0.2885
Residuals                 14 276.0 19.71429               
Total                     21 724.0                        
------------------------------------------------------------------------
CV = 22.2 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.3872083 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------

No significant interaction: analyzing the simple effect
------------------------------------------------------------------------
Variedades
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels Means
1     a1 18.25
2     a2 21.75
------------------------------------------------------------------------
Densidad
T test (LSD)
------------------------------------------------------------------------
Groups  Treatments  Means
a    b2      23.25 
 b   b1      16.75 
------------------------------------------------------------------------

Época
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels Means
1     c1    19
2     c2    21
------------------------------------------------------------------------

DBCA factorial pxqxr con tratamientos adicionales con el paquete ExpDes

data(ex6)
attach(ex6)
data(respAd)
fat3.ad.rbd(factor1 = fatorA,
  factor2 =  fatorB,
  factor3 =  fatorC,
  block =  rep,
  resp =  resp,
  respAd =  respAd,
  quali = c(TRUE,TRUE, TRUE),
  mcomp = "tukey",
  fac.names =
    c("Factor A",
      "Factor B",
      "Factor C"), 
  sigT = 0.05, 
  sigF = 0.05)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Factor A 
FACTOR 2:  Factor B 
FACTOR 3:  Factor C 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                           DF       SS      MS     Fc  Pr>Fc
Block                       2  0.82889 0.41444 0.7397 0.4929
Factor A                    1  0.51042 0.51042  0.911  0.354
Factor B                    1  0.51042 0.51042  0.911  0.354
Factor C                    1  0.05042 0.05042   0.09 0.7681
Factor A*Factor B           1  2.87042 2.87042 5.1232 0.0379
Factor A*Factor C           1  0.12042 0.12042 0.2149 0.6492
Factor B*Factor C           1  0.40042 0.40042 0.7147 0.4104
Factor A*Factor B*Factor C  1  0.02042 0.02042 0.0364  0.851
Ad vs Factorial             1  0.57042 0.57042 1.0181  0.328
Residuals                  16  8.96444 0.56028              
Total                      26 14.84667                      
------------------------------------------------------------------------
CV = 7.39 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.6702706 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
According to the F test, the means of the two groups are statistical equal.
              Means
Additional 10.53333
Factorial  10.07083
------------------------------------------------------------------------



Significant Factor A*Factor B  interaction: analyzing the interaction
------------------------------------------------------------------------

Analyzing  Factor A  inside of each level of  Factor B 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                    DF      SS      MS     Fc  Pr>Fc
Factor A:Factor B 1  1 2.90083 2.90083 5.1775  0.037
Factor A:Factor B 2  1 0.48000 0.48000 0.8567 0.3684
Residuals           16 8.96444 0.56028              
------------------------------------------------------------------------



 Factor A  inside of the level  1  of  Factor B 
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    1   10.41667 
 b   2   9.433333 
------------------------------------------------------------------------


 Factor A  inside of the level  2  of  Factor B 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels    Means
1      1 10.01667
2      2 10.41667
------------------------------------------------------------------------



Analyzing  Factor B  inside of each level of  Factor A 
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                    DF      SS      MS     Fc  Pr>Fc
Factor B:Factor A 1  1 0.48000 0.48000 0.8567 0.3684
Factor B:Factor A 2  1 2.90083 2.90083 5.1775  0.037
Residuals           16 8.96444 0.56028              
------------------------------------------------------------------------



 Factor B  inside of the level  1  of  Factor A 

According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels    Means
1      1 10.41667
2      2 10.01667
------------------------------------------------------------------------


 Factor B  inside of the level  2  of  Factor A 
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a    2   10.41667 
 b   1   9.433333 
------------------------------------------------------------------------

Analizing the effect of the factor  Factor C 
------------------------------------------------------------------------
Factor C
According to the F test, the means of this factor are not different.
------------------------------------------------------------------------
  Niveis   Medias
1      1 10.02500
2      2 10.11667
------------------------------------------------------------------------
fat3.ad.rbd(factor1 = fatorA,
  factor2 =  fatorB,
  factor3 =  fatorC,
  block =  rep,
  resp =  resp,
  respAd =  respAd,
  quali = c(TRUE,TRUE, TRUE),
  mcomp = "tukey",
  fac.names =
    c("Factor A",
      "Factor B",
      "Factor C"), 
  sigT = 0.05, 
  sigF = 0.05,
  unfold = 0)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Factor A 
FACTOR 2:  Factor B 
FACTOR 3:  Factor C 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                           DF       SS      MS     Fc  Pr>Fc
Block                       2  0.82889 0.41444 0.7397 0.4929
Factor A                    1  0.51042 0.51042  0.911  0.354
Factor B                    1  0.51042 0.51042  0.911  0.354
Factor C                    1  0.05042 0.05042   0.09 0.7681
Factor A*Factor B           1  2.87042 2.87042 5.1232 0.0379
Factor A*Factor C           1  0.12042 0.12042 0.2149 0.6492
Factor B*Factor C           1  0.40042 0.40042 0.7147 0.4104
Factor A*Factor B*Factor C  1  0.02042 0.02042 0.0364  0.851
Ad vs Factorial             1  0.57042 0.57042 1.0181  0.328
Residuals                  16  8.96444 0.56028              
Total                      26 14.84667                      
------------------------------------------------------------------------
CV = 7.39 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.6702706 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
According to the F test, the means of the two groups are statistical equal.
              Means
Additional 10.53333
Factorial  10.07083
------------------------------------------------------------------------
attach(ex6)
fat3.ad.rbd(factor1 = fatorA,
  factor2 =  fatorB,
  factor3 =  fatorC,
  block =  rep,
  resp =  resp,
  respAd =  respAd,
  quali = c(TRUE,TRUE, TRUE),
  mcomp = "tukey",
  fac.names =
    c("Factor A",
      "Factor B",
      "Factor C"), 
  sigT = 0.05, 
  sigF = 0.05,
  unfold = 1)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Factor A 
FACTOR 2:  Factor B 
FACTOR 3:  Factor C 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                           DF       SS      MS     Fc  Pr>Fc
Block                       2  0.82889 0.41444 0.7397 0.4929
Factor A                    1  0.51042 0.51042  0.911  0.354
Factor B                    1  0.51042 0.51042  0.911  0.354
Factor C                    1  0.05042 0.05042   0.09 0.7681
Factor A*Factor B           1  2.87042 2.87042 5.1232 0.0379
Factor A*Factor C           1  0.12042 0.12042 0.2149 0.6492
Factor B*Factor C           1  0.40042 0.40042 0.7147 0.4104
Factor A*Factor B*Factor C  1  0.02042 0.02042 0.0364  0.851
Ad vs Factorial             1  0.57042 0.57042 1.0181  0.328
Residuals                  16  8.96444 0.56028              
Total                      26 14.84667                      
------------------------------------------------------------------------
CV = 7.39 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.6702706 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
According to the F test, the means of the two groups are statistical equal.
              Means
Additional 10.53333
Factorial  10.07083
------------------------------------------------------------------------

No significant interaction: analyzing the simple effect
------------------------------------------------------------------------
Factor A
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels    Means
1      1 10.21667
2      2  9.92500
------------------------------------------------------------------------
Factor B
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels    Means
1      1  9.92500
2      2 10.21667
------------------------------------------------------------------------
Factor C
According to the F test, the means of this factor are statistical equal.
------------------------------------------------------------------------
  Levels    Means
1      1 10.02500
2      2 10.11667
------------------------------------------------------------------------
attach(ex6)
fat3.ad.rbd(factor1 = fatorA,
  factor2 =  fatorB,
  factor3 =  fatorC,
  block =  rep,
  resp =  resp,
  respAd =  respAd,
  quali = c(TRUE,TRUE, TRUE),
  mcomp = "tukey",
  fac.names =
    c("Factor A",
      "Factor B",
      "Factor C"), 
  sigT = 0.05, 
  sigF = 0.05,
  unfold = 2)
------------------------------------------------------------------------
Legend:
FACTOR 1:  Factor A 
FACTOR 2:  Factor B 
FACTOR 3:  Factor C 
------------------------------------------------------------------------

------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
                           DF       SS      MS     Fc  Pr>Fc
Block                       2  0.82889 0.41444 0.7397 0.4929
Factor A                    1  0.51042 0.51042  0.911  0.354
Factor B                    1  0.51042 0.51042  0.911  0.354
Factor C                    1  0.05042 0.05042   0.09 0.7681
Factor A*Factor B           1  2.87042 2.87042 5.1232 0.0379
Factor A*Factor C           1  0.12042 0.12042 0.2149 0.6492
Factor B*Factor C           1  0.40042 0.40042 0.7147 0.4104
Factor A*Factor B*Factor C  1  0.02042 0.02042 0.0364  0.851
Ad vs Factorial             1  0.57042 0.57042 1.0181  0.328
Residuals                  16  8.96444 0.56028              
Total                      26 14.84667                      
------------------------------------------------------------------------
CV = 7.39 %

------------------------------------------------------------------------
Shapiro-Wilk normality test
p-value:  0.6702706 
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Contrast of the additional treatment with the factorial
------------------------------------------------------------------------
According to the F test, the means of the two groups are statistical equal.
              Means
Additional 10.53333
Factorial  10.07083
------------------------------------------------------------------------

Prueba de Dunnet

data <- ex6 %>%
  mutate_if(is.integer, factor)%>% 
  dplyr::select(fatorA, fatorB, fatorC, rep, resp) %>%
  dplyr::bind_rows(data.frame(fatorA = "ad",
                       fatorB = "ad",
                       fatorC = "ad",
                       rep = factor(c(1:3)),
                       resp = respAd)) %>%
  dplyr::mutate(tratamiento =
                  interaction(fatorA,
                              fatorB,
                              fatorC))

modelo <- lm(resp ~ tratamiento + rep, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = resp ~ tratamiento + rep, data = data)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)
2.1.1 - 1.1.1 == 0    -0.90000    0.61116  -1.473    0.580
1.2.1 - 1.1.1 == 0    -0.20000    0.61116  -0.327    1.000
2.2.1 - 1.1.1 == 0     0.40000    0.61116   0.654    0.985
1.1.2 - 1.1.1 == 0     0.43333    0.61116   0.709    0.976
2.1.2 - 1.1.1 == 0    -0.63333    0.61116  -1.036    0.858
1.2.2 - 1.1.1 == 0    -0.16667    0.61116  -0.273    1.000
2.2.2 - 1.1.1 == 0     0.03333    0.61116   0.055    1.000
ad.ad.ad - 1.1.1 == 0  0.33333    0.61116   0.545    0.995
(Adjusted p values reported -- single-step method)

Especificar el nivel del tratamiento adicional

data$tratamiento <-
  relevel(data$tratamiento,"ad.ad.ad")
modelo <- lm(resp ~ tratamiento + rep, data = data)
du <- glht(modelo, linfct = mcp(tratamiento = "Dunnet"))
summary(du)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: lm(formula = resp ~ tratamiento + rep, data = data)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)
1.1.1 - ad.ad.ad == 0 -0.33333    0.61116  -0.545    0.995
2.1.1 - ad.ad.ad == 0 -1.23333    0.61116  -2.018    0.274
1.2.1 - ad.ad.ad == 0 -0.53333    0.61116  -0.873    0.932
2.2.1 - ad.ad.ad == 0  0.06667    0.61116   0.109    1.000
1.1.2 - ad.ad.ad == 0  0.10000    0.61116   0.164    1.000
2.1.2 - ad.ad.ad == 0 -0.96667    0.61116  -1.582    0.510
1.2.2 - ad.ad.ad == 0 -0.50000    0.61116  -0.818    0.950
2.2.2 - ad.ad.ad == 0 -0.30000    0.61116  -0.491    0.997
(Adjusted p values reported -- single-step method)

Diseño cuadrado latino en arreglo factorial con tres vías

Planeamiento

Crear un libro de campo con el paquete agricolae

trt <- c(2,2,2)

salida <- agricolae::design.ab(trt = trt,
                               design = "lsd",
                               serie = 3,
                               seed = 123,
                               kinds = "Super-Duper",
                               randomization = TRUE)

salida$book <- salida$book %>%
  dplyr::mutate(A = factor(A,
                           levels = c(1,2),
                           labels = c("Variedad A",
                                      "Variedad B")),
                B = factor(B,
                           levels = c(1,2),
                           labels = c("Dosis 50",
                                      "Dosis 100")),
                C = factor(C,
                           levels = c(1,2),
                           labels = c("Dosis 20",
                                      "Dosis 40"))) %>%
  dplyr::rename("Variedad" = "A",
                "Dosis_N" = "B",
                "Dosis_K" = "C")
salida$book %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)

Guardar el libro generado

write.table(salida$book,
            "books/lsd3w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida$book,
           "books/lsd3w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)

Guardar el libro de campo generado

write.table(salida %>% zigzag(),
            "books/lsd3w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(salida %>% zigzag(),
           "books/lsd3w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
agricolaeplotr::plot_design.factorial_lsd(
  design = salida,
  factor_name = "Variedad",
  reverse_y = TRUE) +
  labs(fill = "Variedades",
       x = "Columnas",
       y = "Filas")
agricolaeplotr::plot_design.factorial_lsd(
  design = salida,
  factor_name = "Dosis_N",
  reverse_y = TRUE) +
  labs(fill = "Dosis de Nitrógeno",
       x = "Columnas",
       y = "Filas")
agricolaeplotr::plot_design.factorial_lsd(
  design = salida,
  factor_name = "Dosis_K",
  reverse_y = TRUE) +
  labs(fill = "Dosis de Potasio",
       x = "Columnas",
       y = "Filas")

Crear un libro de campo con el paquete edibble

menu_factorial()
design("Factorial Design") %>%
  set_units(unit = 280) %>%
  set_trts(trt1 = 4,
           trt2 = 7) %>%
  allot_trts(~unit) %>%
  assign_trts("random", seed = 60) %>%
  serve_table()
lsd <- takeout(menu_factorial(trt = c(2,2,2),
                              r = 8,
                              seed = 123))
lsd %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
lsd2 <- design("Factorial Design") %>%
  set_units(row = 8,
            col = 8,
            unit = crossed_by(row, col)) %>%
  set_trts(trt1 = 2,
           trt2 = 2,
           trt3 = 2) %>%
  allot_trts(trt1 + trt2 + trt3 ~ unit) %>%
  assign_trts("random", seed = 123) %>%
  serve_table()
lsd2 %>% 
  gt::gt() %>%
  gt::opt_interactive(use_search = TRUE,
                      use_filters = TRUE,
                      use_compact_mode = TRUE,
                      page_size_default = 5)
write.table(lsd2 %>% as.data.frame(),
            "books/rcbd3w.txt",
            row.names = FALSE,
            sep = "\t")

write.xlsx(lsd2 %>% as.data.frame(),
           "books/rcbd3w.xlsx",
           sheetName = "book",
           append = FALSE,
           row.names = FALSE)
deggust::autoplot(lsd2)
plot(lsd2)

Análisis de DCL en arreglo factorial con tres vías

Importación de datos

archivos <- list.files(pattern = "datos dcl 3w.xlsx", 
                       full.names = TRUE,
                       recursive = TRUE)

# Importación
data <- readxl::read_xlsx(archivos,
                           sheet = "Hoja1")

# Preprocesamiento

data <- data %>%
  mutate_if(is.character, factor) %>%
  mutate(row = factor(row),
         col = factor(col))

Creación del modelo lineal con interacción

modelo.dcl1 <- lm(rdto ~ A*B*C+row + col, data = data)
modelo.dcl2 <- lm(rdto ~ A*B+C+row + col, data = data)
modelo.dcl3 <- lm(rdto ~ A*B+C+row + col, data = data)
modelo.dcl4 <- lm(rdto ~ A+B*C+row + col, data = data)
modelo.dcl5 <- lm(rdto ~ A+B+C+row + col, data = data)
modelo.dcl6 <- lm(rdto ~ A*B*C+row, data = data)
modelo.dcl7 <- lm(rdto ~ A*B+C+row, data = data)
modelo.dcl8 <- lm(rdto ~ A*B+C+row, data = data)
modelo.dcl9 <- lm(rdto ~ A+B*C+row, data = data)
modelo.dcl10 <- lm(rdto ~ A+B+C+row, data = data)
modelo.dcl11 <- lm(rdto ~ A*B*C+col, data = data)
modelo.dcl12 <- lm(rdto ~ A*B+C+col, data = data)
modelo.dcl13 <- lm(rdto ~ A*B+C+col, data = data)
modelo.dcl14 <- lm(rdto ~ A+B*C+col, data = data)
modelo.dcl15 <- lm(rdto ~ A+B+C+col, data = data)
broom::glance(modelo.dcl1) %>%
  bind_rows(broom::glance(modelo.dcl2),
            broom::glance(modelo.dcl3),
            broom::glance(modelo.dcl4),
            broom::glance(modelo.dcl5),
            broom::glance(modelo.dcl6),
            broom::glance(modelo.dcl7),
            broom::glance(modelo.dcl8),
            broom::glance(modelo.dcl9),
            broom::glance(modelo.dcl10),
            broom::glance(modelo.dcl11),
            broom::glance(modelo.dcl12),
            broom::glance(modelo.dcl13),
            broom::glance(modelo.dcl14),
            broom::glance(modelo.dcl15)) %>%
  dplyr::mutate(Modelo = c("a*b*c+row+col", "a*b+c+row+col",
                           "a*c+b+row+col",
                           "a+b*c+row+col", "a+b+c+row+col",
                           "a*b*c+row", "a*b+c+row",
                           "a*c+b+row",
                           "a+b*c+row", "a+b+c+row",
                           "a*b*c+col", "a*b+c+col",
                           "a*c+b+col",
                           "a+b*c+col", "a+b+c+col")) %>%
  dplyr::select(Modelo, AIC, BIC) %>%
  dplyr::arrange(BIC) %>%
  dplyr::mutate(Mérito = 1:n()) %>%
  dplyr::relocate(Mérito, Modelo) %>%
  gt()
Mérito Modelo AIC BIC
1 a*b+c+row 129.3267 157.3922
2 a*c+b+row 129.3267 157.3922
3 a*b+c+col 130.5685 158.6340
4 a*c+b+col 130.5685 158.6340
5 a*b*c+row 124.4155 158.9576
6 a*b*c+col 125.8854 160.4275
7 a+b+c+row 134.7738 160.6804
8 a+b+c+col 135.8803 161.7869
9 a+b*c+row 136.3099 164.3754
10 a+b*c+col 137.4244 165.4899
11 a*b+c+row+col 137.4279 180.6055
12 a*c+b+row+col 137.4279 180.6055
13 a*b*c+row+col 131.3575 181.0118
14 a+b+c+row+col 143.5501 184.5689
15 a+b*c+row+col 145.0466 188.2243
 
modelo.dcl <- lm(rdto ~ A * B + C + row, data = data)

Definición del modelo

Con filas y columnas

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2 + \beta_4*F_{2} + \beta_5*F_{3} + \beta_6*F_{4} + \beta_7*F_{5} + \beta_8*F_{6} + \beta_9*F_{7} + \beta_{10}*F_{8} + \beta_{11}*C_2 + \beta_{12}*C_3 + \beta_{13}*C_4 + \beta_{14}*C_5 \beta_{15}*C_6 \beta_{16}*C_8 \beta_{17}*C_9 + \beta_{18}*A_2*B_2 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2 + \beta_4*F_{2} + \beta_5*F_{3} + \beta_6*F_{4} + \beta_7*F_{5} + \beta_8*F_{6} + \beta_9*F_{7} + \beta_{10}*F_{8} + \beta_{11}*C_2 + \beta_{12}*C_3 + \beta_{13}*C_4 + \beta_{14}*C_5 \beta_{15}*C_6 \beta_{16}*C_8 \beta_{17}*C_9 + \beta_{18}*A_2*B_2\]

Con solo filas

\[Y_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2 + \beta_4*F_{2} + \beta_5*F_{3} + \beta_6*F_{4} + \beta_7*F_{5} + \beta_8*F_{6} + \beta_9*F_{7} + \beta_{10}*F_{8} + \beta_{11}*A_2*B_2 + \epsilon_i\]

\[\hat{Y}_i = \beta_0 + \beta_1*A_2 + \beta_2*B_2 + \beta_3*C_2 + \beta_4*F_{2} + \beta_5*F_{3} + \beta_6*F_{4} + \beta_7*F_{5} + \beta_8*F_{6} + \beta_9*F_{7} + \beta_{10}*F_{8} + \beta_{11}*A_2*B_2\]

summary(modelo.dcl)

Call:
lm(formula = rdto ~ A * B + C + row, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4352 -0.3981  0.0327  0.4434  1.0245 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.01525    0.26057  19.247  < 2e-16 ***
AA2          0.98113    0.21276   4.611 2.64e-05 ***
BB2         -0.44806    0.21276  -2.106   0.0401 *  
CC2          0.38159    0.15044   2.536   0.0142 *  
row2         0.21150    0.30089   0.703   0.4852    
row3         0.25900    0.30089   0.861   0.3933    
row4         0.25462    0.30089   0.846   0.4013    
row5         0.24813    0.30089   0.825   0.4133    
row6         0.29600    0.30089   0.984   0.3298    
row7        -0.00375    0.30089  -0.012   0.9901    
row8        -0.24188    0.30089  -0.804   0.4251    
AA2:BB2      0.76219    0.30089   2.533   0.0144 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6018 on 52 degrees of freedom
Multiple R-squared:  0.6593,    Adjusted R-squared:  0.5872 
F-statistic: 9.149 on 11 and 52 DF,  p-value: 8.098e-09

Verificación visual de los supuestos del modelo

performance::check_model(modelo.dcl)
ggResidpanel::resid_panel(modelo.dcl)
influence.measures(modelo.dcl)
Influence measures of
     lm(formula = rdto ~ A * B + C + row, data = data) :

      dfb.1_   dfb.AA2   dfb.BB2   dfb.CC2  dfb.row2  dfb.row3  dfb.row4
1   2.66e-01 -1.30e-01 -1.30e-01  0.092088 -1.84e-01 -1.84e-01 -1.84e-01
2   3.25e-01  1.03e-15  9.79e-16 -0.140519 -2.81e-01 -2.81e-01 -2.81e-01
3   4.00e-01 -1.63e-01 -1.63e-01 -0.115405 -2.31e-01 -2.31e-01 -2.31e-01
4   8.78e-03  7.17e-03 -3.68e-18  0.005072 -1.01e-02 -1.01e-02 -1.01e-02
5   7.24e-02  4.80e-17  2.46e-18  0.041789 -8.36e-02 -8.36e-02 -8.36e-02
6  -3.27e-01 -2.01e-01  1.59e-16  0.141791  2.84e-01  2.84e-01  2.84e-01
7  -9.02e-01 -9.58e-17 -5.52e-01  0.390612  7.81e-01  7.81e-01  7.81e-01
8   1.94e-01 -1.21e-17  1.59e-01  0.112170 -2.24e-01 -2.24e-01 -2.24e-01
9  -1.03e-02  1.12e-17  3.15e-18  0.017829  3.57e-02  1.58e-18  2.55e-18
10  1.79e-17 -7.89e-02  3.83e-17  0.055757 -1.12e-01 -5.92e-17 -5.81e-17
11  5.84e-17  1.11e-17 -1.30e-17 -0.094614  1.89e-01 -6.51e-17 -6.03e-17
12  1.27e-01  1.89e-16 -3.10e-01 -0.219192 -4.38e-01 -1.42e-16 -1.53e-16
13  8.26e-03 -2.02e-02 -2.81e-19 -0.014299 -2.86e-02 -2.80e-18 -2.91e-18
14  1.81e-17  1.17e-16 -1.45e-01  0.102757 -2.06e-01 -5.86e-17 -6.88e-17
15  1.66e-01 -2.03e-01 -2.03e-01 -0.143586  2.87e-01 -6.69e-17 -6.08e-17
16  7.69e-02 -1.88e-01 -1.88e-01  0.133266  2.67e-01 -1.74e-17 -1.54e-17
17  9.94e-03 -2.43e-02 -4.22e-18 -0.017211 -1.23e-17 -3.44e-02 -6.00e-18
18  4.21e-18  3.09e-17 -5.31e-02  0.037546 -2.52e-17 -7.51e-02 -1.03e-17
19 -8.25e-18  9.10e-02 -9.47e-18 -0.064350  5.13e-17  1.29e-01  2.90e-17
20 -5.10e-02  1.25e-01  1.25e-01 -0.088419 -8.90e-17 -1.77e-01 -3.90e-17
21 -1.02e-01 -9.73e-17  2.50e-01  0.177008  1.41e-16  3.54e-01  5.35e-17
22 -2.55e-01  3.13e-01  3.13e-01  0.221089 -1.58e-16 -4.42e-01 -6.29e-17
23  1.07e-20  2.43e-19  8.70e-21 -0.000886  1.12e-19  1.77e-03 -3.55e-19
24 -6.93e-02  8.49e-17  2.95e-17  0.120100  2.52e-17  2.40e-01 -2.85e-17
25  1.23e-17  6.78e-02 -3.34e-17 -0.047908  2.10e-17  2.03e-17  9.58e-02
26  9.52e-02 -2.33e-01 -2.33e-01  0.164920  4.65e-17  3.50e-17  3.30e-01
27 -1.55e-01 -2.69e-16  3.80e-01  0.268516  1.63e-16  1.14e-16  5.37e-01
28  4.42e-03 -5.42e-03 -5.42e-03 -0.003829  1.98e-18  1.46e-18  7.66e-03
29 -5.45e-17  1.86e-16 -2.35e-01  0.166050 -9.27e-17 -6.34e-17 -3.32e-01
30  1.56e-01  1.59e-17  1.25e-16 -0.270524  9.26e-17  9.18e-17 -5.41e-01
31 -1.67e-02  4.10e-02 -4.27e-18  0.028978  6.96e-18  6.14e-18  5.80e-02
32 -7.03e-17  2.60e-17  4.21e-17  0.078000  3.56e-17  4.30e-17 -1.56e-01
33 -1.21e-16  1.83e-16 -2.03e-01  0.143341 -4.45e-18  2.84e-17  2.47e-17
34  1.49e-02 -6.61e-18  4.32e-18 -0.025889  2.35e-17  2.71e-17  3.12e-17
35 -1.11e-01  2.73e-01  2.73e-01 -0.192935 -2.18e-17  2.47e-17  4.98e-17
36  2.62e-16 -3.32e-17 -6.10e-17 -0.188231 -1.72e-16 -2.05e-16 -2.27e-16
37 -2.25e-01  2.75e-01  2.75e-01  0.194651  1.98e-17  5.52e-17  8.37e-17
38  4.37e-02 -1.07e-01  6.69e-18 -0.075693  4.63e-18  1.29e-17  1.95e-17
39 -1.05e-01 -1.40e-16  2.58e-01  0.182291  3.04e-18 -3.88e-17 -4.70e-17
40  2.76e-16  3.80e-01 -1.71e-16 -0.268771 -6.35e-17 -6.87e-17 -1.39e-16
41 -7.32e-17  5.78e-18  3.47e-17  0.147311  2.16e-17  5.37e-17  5.84e-17
42 -9.20e-02 -1.25e-16  2.25e-01  0.159306  9.70e-17  4.57e-17  5.76e-17
43  4.67e-03 -1.14e-02 -6.35e-19 -0.008094 -4.23e-18 -2.89e-18 -2.23e-18
44  5.92e-17 -1.70e-16  2.27e-01 -0.160514  5.27e-17  1.14e-17  2.76e-17
45  5.04e-17  1.85e-01 -8.22e-17 -0.130918  4.54e-17  2.60e-17  2.25e-17
46  3.69e-02 -9.05e-02 -9.05e-02  0.063975  3.17e-17  1.65e-17  1.21e-17
47  1.83e-01 -7.46e-17  4.98e-17 -0.316897  7.08e-17  1.27e-16  1.37e-16
48 -6.21e-02  7.61e-02  7.61e-02  0.053818 -2.41e-17 -1.16e-17 -1.02e-17
49  4.64e-02 -5.68e-02 -5.68e-02 -0.040163 -6.92e-18 -1.57e-17 -1.41e-17
50  5.82e-02 -1.43e-01  9.90e-18 -0.100836  5.52e-18  1.71e-17  1.58e-17
51  6.91e-02 -4.70e-18  2.11e-17 -0.119696  8.54e-17  9.00e-17  1.01e-16
52  1.10e-17  1.62e-02 -1.01e-17 -0.011489 -1.80e-18 -2.63e-18 -3.30e-18
53  1.57e-16 -3.83e-17 -4.02e-17 -0.097503 -1.16e-16 -1.28e-16 -1.32e-16
54  7.82e-02  1.22e-16 -1.92e-01 -0.135446 -7.91e-17 -4.27e-17 -6.43e-17
55  3.91e-02 -9.58e-02 -9.58e-02  0.067735  2.82e-18 -1.06e-17 -9.74e-18
56 -1.33e-17 -2.07e-16  1.96e-01 -0.138449  1.13e-16  7.97e-17  9.99e-17
57 -6.25e-02 -9.03e-17  1.53e-01  0.108269  1.08e-16  7.88e-17  9.92e-17
58  3.98e-03 -4.88e-03 -4.88e-03 -0.003450 -6.39e-19 -1.54e-18 -1.15e-18
59  3.66e-17  1.73e-16 -1.82e-01  0.128759 -1.12e-16 -7.69e-17 -1.07e-16
60 -4.97e-02  1.61e-17 -5.91e-18  0.086098 -3.34e-17 -4.28e-17 -4.41e-17
61 -4.87e-02  1.19e-01  1.19e-01 -0.084354 -1.57e-17  1.56e-18 -2.01e-18
62  1.70e-16 -3.73e-17 -3.43e-17 -0.151993 -1.19e-16 -1.45e-16 -1.37e-16
63  9.65e-18  5.32e-02 -3.21e-17 -0.037625  2.02e-17  1.64e-17  1.70e-17
64  1.01e-01 -2.48e-01 -2.93e-17 -0.175649 -7.14e-17 -4.73e-17 -5.76e-17
    dfb.row5  dfb.row6  dfb.row7  dfb.row8  dfb.AA2.    dffit cov.r   cook.d
1  -1.84e-01 -1.84e-01 -1.84e-01 -1.84e-01  0.092088  0.31900 1.401 8.57e-03
2  -2.81e-01 -2.81e-01 -2.81e-01 -2.81e-01  0.140519  0.48677 1.223 1.97e-02
3  -2.31e-01 -2.31e-01 -2.31e-01 -2.31e-01  0.115405  0.39978 1.322 1.34e-02
4  -1.01e-02 -1.01e-02 -1.01e-02 -1.01e-02 -0.005072  0.01757 1.553 2.62e-05
5  -8.36e-02 -8.36e-02 -8.36e-02 -8.36e-02  0.041789  0.14476 1.521 1.78e-03
6   2.84e-01  2.84e-01  2.84e-01  2.84e-01  0.141791 -0.49118 1.218 2.01e-02
7   7.81e-01  7.81e-01  7.81e-01  7.81e-01  0.390612 -1.35312 0.274 1.35e-01
8  -2.24e-01 -2.24e-01 -2.24e-01 -2.24e-01 -0.112170  0.38857 1.333 1.27e-02
9   1.19e-18  6.73e-18  7.00e-18 -7.42e-18  0.017829  0.06176 1.548 3.24e-04
10 -4.60e-17 -5.20e-17 -5.19e-17 -8.51e-17  0.055757 -0.19315 1.496 3.16e-03
11 -7.35e-17 -8.63e-17 -7.87e-17 -7.09e-17  0.094614  0.32775 1.393 9.04e-03
12 -1.09e-16 -1.26e-16 -1.13e-16 -1.10e-16  0.219192 -0.75930 0.875 4.67e-02
13 -6.85e-19 -9.28e-19 -2.21e-18  3.97e-19  0.014299 -0.04953 1.550 2.08e-04
14 -4.38e-17 -6.16e-17 -3.58e-17 -3.71e-17  0.102757 -0.35596 1.366 1.07e-02
15 -7.42e-17 -7.42e-17 -7.11e-17 -1.08e-16  0.143586  0.49740 1.210 2.06e-02
16 -2.78e-17 -3.05e-17 -2.08e-17 -4.81e-17  0.133266  0.46165 1.253 1.78e-02
17 -6.88e-18 -8.76e-18 -1.05e-17 -5.25e-18  0.017211 -0.05962 1.548 3.02e-04
18 -1.10e-17 -1.64e-17 -1.07e-17 -5.21e-18  0.037546 -0.13006 1.527 1.44e-03
19  3.07e-17  3.63e-17  4.28e-17  2.68e-17 -0.064350  0.22291 1.477 4.20e-03
20 -5.55e-17 -5.32e-17 -6.79e-17 -4.66e-17 -0.088419 -0.30629 1.413 7.91e-03
21  5.63e-17  8.35e-17  5.93e-17  5.90e-17 -0.177008  0.61318 1.065 3.10e-02
22 -8.48e-17 -1.04e-16 -1.13e-16 -5.52e-17 -0.221089 -0.76587 0.867 4.75e-02
23 -2.70e-19 -3.87e-19 -2.48e-19 -3.69e-19  0.000886  0.00307 1.554 8.01e-07
24 -1.28e-17 -2.98e-17 -8.78e-18 -2.33e-17  0.120100  0.41604 1.304 1.45e-02
25  9.26e-18  1.54e-17  1.54e-17  1.46e-17 -0.047908  0.16596 1.511 2.33e-03
26 -4.30e-18 -4.48e-19  2.21e-17  4.58e-18  0.164920  0.57130 1.119 2.70e-02
27  1.04e-16  1.39e-16  9.39e-17  1.49e-16 -0.268516  0.93017 0.663 6.85e-02
28  6.82e-19  9.52e-19  1.20e-18  7.44e-19  0.003829  0.01327 1.553 1.50e-05
29 -3.45e-17 -7.55e-17 -3.52e-17 -7.83e-17  0.166050 -0.57521 1.114 2.73e-02
30  1.61e-16  1.88e-16  1.61e-16  1.28e-16 -0.270524 -0.93712 0.655 6.94e-02
31 -1.13e-18 -6.71e-20  3.30e-18  2.41e-18 -0.028978  0.10038 1.538 8.55e-04
32  5.49e-17  6.15e-17  6.05e-17  4.98e-17 -0.078000 -0.27020 1.443 6.17e-03
33 -2.87e-01 -1.57e-17  1.52e-17 -3.98e-18  0.143341 -0.49655 1.211 2.05e-02
34 -5.18e-02  2.56e-17  2.40e-17  1.94e-17 -0.025889 -0.08968 1.541 6.83e-04
35 -3.86e-01  2.28e-17  8.75e-18  2.68e-17 -0.192935 -0.66835 0.994 3.66e-02
36  3.76e-01 -1.78e-16 -1.67e-16 -1.62e-16  0.188231  0.65205 1.015 3.49e-02
37 -3.89e-01  4.09e-17  4.77e-17  4.86e-17 -0.194651 -0.67429 0.986 3.72e-02
38 -1.51e-01  4.14e-18  7.07e-18  1.05e-17  0.075693 -0.26221 1.449 5.81e-03
39  3.65e-01 -9.94e-18 -4.45e-17 -1.01e-17 -0.182291  0.63147 1.042 3.28e-02
40  5.38e-01 -9.24e-17 -9.57e-17 -1.27e-16 -0.268771  0.93105 0.662 6.86e-02
41  3.88e-17 -2.95e-01  7.56e-17  4.50e-17 -0.147311 -0.51030 1.195 2.16e-02
42  1.40e-17  3.19e-01  4.36e-17  7.96e-17 -0.159306  0.55185 1.144 2.52e-02
43 -7.10e-19 -1.62e-02 -3.50e-18 -4.27e-18  0.008094 -0.02804 1.552 6.68e-05
44  0.00e+00  3.21e-01  9.67e-18  5.79e-17 -0.160514  0.55604 1.138 2.56e-02
45  0.00e+00  2.62e-01  4.20e-17  4.00e-17 -0.130918  0.45351 1.262 1.72e-02
46  5.62e-18  1.28e-01  2.10e-17  2.66e-17  0.063975  0.22162 1.478 4.16e-03
47  1.39e-16 -6.34e-01  1.30e-16  6.16e-17 -0.316897 -1.09776 0.482 9.29e-02
48 -4.72e-18 -1.08e-01 -2.03e-17 -2.24e-17 -0.053818 -0.18643 1.500 2.94e-03
49 -1.50e-17 -3.64e-18  8.03e-02 -3.34e-18  0.040163  0.13913 1.523 1.64e-03
50  3.54e-17  0.00e+00 -2.02e-01 -3.08e-17  0.100836 -0.34930 1.373 1.03e-02
51  9.80e-17  9.77e-17 -2.39e-01  6.31e-17 -0.119696 -0.41464 1.306 1.44e-02
52 -4.55e-18 -2.08e-18  2.30e-02 -3.19e-19 -0.011489  0.03980 1.551 1.35e-04
53 -1.17e-16 -1.24e-16  1.95e-01 -1.06e-16  0.097503  0.33776 1.384 9.60e-03
54 -2.72e-17 -6.14e-17 -2.71e-01 -1.13e-16  0.135446 -0.46920 1.244 1.84e-02
55 -1.34e-17 -6.14e-18  1.35e-01 -9.40e-18  0.067735  0.23464 1.469 4.66e-03
56  8.90e-17  1.13e-16  2.77e-01  8.07e-17 -0.138449  0.47960 1.232 1.92e-02
57  7.55e-17  8.38e-17  1.04e-16  2.17e-01 -0.108269  0.37505 1.347 1.18e-02
58 -1.29e-18 -8.47e-19 -6.63e-19  6.90e-03  0.003450  0.01195 1.554 1.21e-05
59 -8.75e-17 -9.48e-17 -7.43e-17 -2.58e-01  0.128759 -0.44604 1.271 1.66e-02
60 -4.14e-17 -5.56e-17 -4.97e-17  1.72e-01  0.086098  0.29825 1.419 7.50e-03
61  3.91e-18  2.24e-18 -1.62e-17 -1.69e-01 -0.084354 -0.29221 1.425 7.20e-03
62 -1.16e-16 -1.59e-16 -1.46e-16  3.04e-01  0.151993  0.52652 1.175 2.30e-02
63  1.05e-17  1.31e-17  2.89e-17  7.53e-02 -0.037625  0.13034 1.527 1.44e-03
64 -1.51e-17 -5.44e-17 -6.76e-17 -3.51e-01  0.175649 -0.60846 1.071 3.05e-02
     hat inf
1  0.188    
2  0.187    
3  0.188    
4  0.188    
5  0.187    
6  0.187    
7  0.188   *
8  0.188    
9  0.188    
10 0.188    
11 0.188    
12 0.188    
13 0.188    
14 0.188    
15 0.188    
16 0.188    
17 0.188    
18 0.187    
19 0.188    
20 0.188    
21 0.187    
22 0.188    
23 0.188    
24 0.188    
25 0.188    
26 0.187    
27 0.187    
28 0.188    
29 0.187    
30 0.187    
31 0.188    
32 0.187    
33 0.187    
34 0.187    
35 0.187    
36 0.187    
37 0.187    
38 0.187    
39 0.187    
40 0.187    
41 0.187    
42 0.187    
43 0.188    
44 0.187    
45 0.188    
46 0.187    
47 0.187    
48 0.188    
49 0.187    
50 0.187    
51 0.187    
52 0.187    
53 0.187    
54 0.188    
55 0.188    
56 0.188    
57 0.188    
58 0.187    
59 0.188    
60 0.188    
61 0.187    
62 0.187    
63 0.187    
64 0.188    
influenceIndexPlot(modelo.dcl)

Cumplimiento de supuestos del modelo lineal general

Independencia de residuos

\(H_0: \text{Los residuos del rendimiento son completamente aleatorios e independientes}\)

\(H_1: \text{Los residuos del rendimiento no son completamente aleatorios e independientes}\)

durbinWatsonTest(modelo.dcl,
                 reps = 5000,
                 max.lag = 10)
 lag Autocorrelation D-W Statistic p-value
   1     -0.01688919      2.002034  0.3892
   2     -0.04870879      2.048467  0.6892
   3     -0.02901524      1.979496  0.6400
   4     -0.05379072      2.023173  0.9248
   5     -0.05581067      2.019674  0.7644
   6     -0.03552560      1.949277  0.8208
   7     -0.11383880      1.996510  0.5080
   8     -0.06219868      1.873342  0.6180
   9     -0.20501407      2.143135  0.1112
  10     -0.02565903      1.778074  0.6448
 Alternative hypothesis: rho[lag] != 0
dwtest(modelo.dcl, alternative = "two.sided")

    Durbin-Watson test

data:  modelo.dcl
DW = 2.002, p-value = 0.3869
alternative hypothesis: true autocorrelation is not 0

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.

Normalidad de residuos

\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)

shapiro.test(rstudent(modelo.dcl))

    Shapiro-Wilk normality test

data:  rstudent(modelo.dcl)
W = 0.97724, p-value = 0.2833
ad.test(rstudent(modelo.dcl))

    Anderson-Darling normality test

data:  rstudent(modelo.dcl)
A = 0.44959, p-value = 0.2683
lillie.test(rstudent(modelo.dcl))

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  rstudent(modelo.dcl)
D = 0.072802, p-value = 0.5461
ks.test(rstudent(modelo.dcl), "pnorm",
        alternative = "two.sided")

    Exact one-sample Kolmogorov-Smirnov test

data:  rstudent(modelo.dcl)
D = 0.076405, p-value = 0.8214
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dcl))

    Cramer-von Mises normality test

data:  rstudent(modelo.dcl)
W = 0.070988, p-value = 0.267
pearson.test(rstudent(modelo.dcl))

    Pearson chi-square normality test

data:  rstudent(modelo.dcl)
P = 13.344, p-value = 0.1006
sf.test(rstudent(modelo.dcl))

    Shapiro-Francia normality test

data:  rstudent(modelo.dcl)
W = 0.97868, p-value = 0.282

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.

Homocedasticidad

\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento

\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento

ncvTest(modelo.dcl)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.08116, Df = 1, p = 0.29844
bptest(modelo.dcl)

    studentized Breusch-Pagan test

data:  modelo.dcl
BP = 12.585, df = 11, p-value = 0.3213
bptest(modelo.dcl, studentize = F)

    Breusch-Pagan test

data:  modelo.dcl
BP = 9.9237, df = 11, p-value = 0.5373
olsrr::ols_test_breusch_pagan(modelo.dcl)

 Breusch Pagan Test for Heteroskedasticity
 -----------------------------------------
 Ho: the variance is constant            
 Ha: the variance is not constant        

              Data               
 --------------------------------
 Response : rdto 
 Variables: fitted values of rdto 

        Test Summary         
 ----------------------------
 DF            =    1 
 Chi2          =    1.08116 
 Prob > Chi2   =    0.2984382

Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.

Recomendación. Debido a que no se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.

Estadísticas globales

modelo.dcl %>% gvlma()

Call:
lm(formula = rdto ~ A * B + C + row, data = data)

Coefficients:
(Intercept)          AA2          BB2          CC2         row2         row3  
    5.01525      0.98113     -0.44806      0.38159      0.21150      0.25900  
       row4         row5         row6         row7         row8      AA2:BB2  
    0.25463      0.24813      0.29600     -0.00375     -0.24187      0.76219  


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = .) 

                     Value  p-value                   Decision
Global Stat        10.6436 0.030875 Assumptions NOT satisfied!
Skewness            1.6596 0.197661    Assumptions acceptable.
Kurtosis            0.4771 0.489754    Assumptions acceptable.
Link Function       8.3284 0.003903 Assumptions NOT satisfied!
Heteroscedasticity  0.1786 0.672568    Assumptions acceptable.

Análisis de varianza

\[Y_{ijklm} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_{l} + \omega_{m} + \epsilon_{ijklm}\]

\[\hat{Y}_{ijklm} = \mu + \tau_{i} + \beta_{j} + \gamma_{k} + \delta_{l} + \omega_{m}\]

Dónde:

\(Y_{ijklm}\) = Valor observado de la variable respuesta.

\(\hat{Y}_{ijklm}\) = Valor ajustado de la variable respuesta.

\(\mu\) = Promedio observado de la variable respuesta.

\(\tau_{i}\) = Efecto del i-ésimo nivel del factor A.

\(\beta_{j}\) = Efecto del j-ésimo nivel del factor B.

\(\gamma_{k}\) = Efecto del k-ésimo nivel del factor C.

\(\delta_{l}\) = Efecto del l-ésimo nivel del factor Fila.

\(\omega_{m}\) = Efecto del m-ésimo nivel del factor Columna.

\(\epsilon_{ijklm}\) = Residuo observado del modelo.

Pruebas de hipótesis

Para el factor A:

\(H_0: \tau_{A1} = \tau_{A2} = 0\)

\(H_1: \text{En al menos un nivel del factor A el } \tau \text{ es diferente a los demás.}\)

\(H_1: \tau_i \neq 0\text{; en al menos un nivel del factor A.}\)

Para el factor B:

\(H_0: \beta_{B1} = \beta_{B2} = 0\)

\(H_1: \text{En al menos un nivel del factor B el } \beta \text{ es diferente a los demás.}\)

\(H_1: \beta_{j} \neq 0\text{; en al menos un nivel del factor B.}\)

Para el factor C:

\(H_0: \gamma_{C1} = \gamma_{C2} = 0\)

\(H_1: \text{En al menos un nivel del factor C el } \gamma \text{ es diferente a los demás.}\)

\(H_1: \gamma_{k} \neq 0\text{; en al menos un nivel del factor C.}\)

Para la interacción entre factor A y factor B:

\(H_0: (\tau\beta)_{A1B1} = (\tau\beta)_{A1B2} = (\tau\beta)_{A2B1} = (\tau\beta)_{A2B2} = 0\)

\(H_1: \text{En al menos una interacción entre el factor A y el factor B el } (\tau\beta) \text{ es diferente a los demás.}\)

\(H_1: (\tau\beta)_{ij} \neq 0\text{; en al menos una interacción entre el factor A y el factor B.}\)

anova(modelo.dcl, test = "F")
Analysis of Variance Table

Response: rdto
          Df  Sum Sq Mean Sq F value    Pr(>F)    
A          1 29.6902 29.6902 81.9882 2.857e-12 ***
B          1  0.0718  0.0718  0.1982   0.65806    
C          1  2.3298  2.3298  6.4337   0.01424 *  
row        7  2.0270  0.2896  0.7996   0.59123    
A:B        1  2.3237  2.3237  6.4168   0.01436 *  
Residuals 52 18.8307  0.3621                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Valor de la tabla de F para el factor A con una significancia de 0.05.

qf(0.95, 1, 45)
[1] 4.056612

Valor de la tabla de F para el factor B con una significancia de 0.05.

qf(0.95, 1, 45)
[1] 4.056612

Valor de la tabla de F para el factor C con una significancia de 0.05.

qf(0.95, 1, 45)
[1] 4.056612

Valor de la tabla de F para la interacción AxB con una significancia de 0.05.

qf(0.95, 1, 45)
[1] 4.056612

Conclusión.

Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor A tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor A.

Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor B tienen un efecto sobre el rendimiento estadísticamente similar a 0.

Con respecto al Factor C: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor C tiene un efecto sobre el rendimiento estadísticamente diferente del resto de niveles del factor C.

agricolae::cv.model(modelo.dcl)
[1] 10.39134

Comparaciones de medias

Para los niveles del factor A:

  • A1 vs A2:

\(H_0: \mu_{A1} - \mu_{A2} = 0\)

\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)

Para los niveles del factor C:

  • C1 vs C2:

\(H_0: \mu_{C1} - \mu_{C2} = 0\)

\(H_1: \mu_{C1} - \mu_{C2} \neq 0\)

Prueba de LSD

Para los niveles del Factor A:

agricolae::LSD.test(modelo.dcl, trt = "A", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "A"

LSD t Test for rdto 

Mean Square Error:  0.3621281 

A,  means and individual ( 95 %) CI

       rdto       std  r      LCL      UCL   Min   Max
A1 5.109969 0.7599529 32 4.896504 5.323434 3.132 6.294
A2 6.472188 0.4977219 32 6.258722 6.685653 5.442 7.422

Alpha: 0.05 ; DF Error: 52
Critical Value of t: 2.006647 

least Significant Difference: 0.3018854 

Treatments with the same letter are not significantly different.

       rdto groups
A2 6.472188      a
A1 5.109969      b
agricolae::LSD.test(modelo.dcl, trt = "A", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "A"

LSD t Test for rdto 

Mean Square Error:  0.3621281 

A,  means and individual ( 95 %) CI

       rdto       std  r      LCL      UCL   Min   Max
A1 5.109969 0.7599529 32 4.896504 5.323434 3.132 6.294
A2 6.472188 0.4977219 32 6.258722 6.685653 5.442 7.422

Alpha: 0.05 ; DF Error: 52
Critical Value of t: 2.006647 

Comparison between treatments means

        difference pvalue signif.       LCL       UCL
A1 - A2  -1.362219      0     *** -1.664104 -1.060333

Para los niveles del Factor C:

agricolae::LSD.test(modelo.dcl, trt = "C", alpha = 0.05,
         group = TRUE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "C"

LSD t Test for rdto 

Mean Square Error:  0.3621281 

C,  means and individual ( 95 %) CI

       rdto       std  r      LCL      UCL   Min   Max
C1 5.600281 1.0739613 32 5.386816 5.813746 3.132 7.289
C2 5.981875 0.7446186 32 5.768410 6.195340 4.315 7.422

Alpha: 0.05 ; DF Error: 52
Critical Value of t: 2.006647 

least Significant Difference: 0.3018854 

Treatments with the same letter are not significantly different.

       rdto groups
C2 5.981875      a
C1 5.600281      b
agricolae::LSD.test(modelo.dcl, trt = "C", alpha = 0.05,
         group = FALSE, main = NULL, console = TRUE)

Study: modelo.dcl ~ "C"

LSD t Test for rdto 

Mean Square Error:  0.3621281 

C,  means and individual ( 95 %) CI

       rdto       std  r      LCL      UCL   Min   Max
C1 5.600281 1.0739613 32 5.386816 5.813746 3.132 7.289
C2 5.981875 0.7446186 32 5.768410 6.195340 4.315 7.422

Alpha: 0.05 ; DF Error: 52
Critical Value of t: 2.006647 

Comparison between treatments means

        difference pvalue signif.        LCL         UCL
C1 - C2 -0.3815938 0.0142       * -0.6834791 -0.07970838

Comparaciones de medias para las interacciones

Para los niveles del factor A dentro del nivel B1:

  • A1 vs A2:

\(H_0: \mu_{A1} - \mu_{A2} = 0\)

\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)

Para los niveles del factor A dentro del nivel B2:

  • A1 vs A2:

\(H_0: \mu_{A1} - \mu_{A2} = 0\)

\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)

Para los niveles del factor B dentro del nivel A1:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

Para los niveles del factor B dentro del nivel A2:

  • B1 vs B2:

\(H_0: \mu_{B1} - \mu_{B2} = 0\)

\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)

Análisis de varianza para interacción de dos factores con el paquete phia

Comparación de los niveles de B dentro de cada nivel de A

phia::testInteractions(modelo.dcl,
                       fixed = "A",
                       across = "B",
                       adjustment = "none")
F Test: 
P-value adjustment method: none
             Value Df Sum of Sq      F  Pr(>F)  
A1         0.44806  1    1.6061 4.4351 0.04005 *
A2        -0.31412  1    0.7894 2.1799 0.14586  
Residuals          52   18.8307                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparación de los niveles de A dentro de cada nivel de B

phia::testInteractions(modelo.dcl,
                       fixed = "B",
                       across = "A",
                       adjustment = "none")
F Test: 
P-value adjustment method: none
             Value Df Sum of Sq      F    Pr(>F)    
B1        -0.98113  1    7.7009 21.265 2.638e-05 ***
B2        -1.74331  1   24.3131 67.139 6.295e-11 ***
Residuals          52   18.8307                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot de interacciones

phia::interactionMeans(model = modelo.dcl,
                       factors = c("A","B")) %>%
  plot()

Comparaciones de medias de los niveles de B dentro de cada nivel de A

FA_A1 <- data %>% filter(A=="A1")
FA_A2 <- data %>% filter(A=="A2")

fcomp1 <- function(x){
  comp <- LSD.test(x$rdto, # Cambiar según nombre de variable respuesta
                   x$B, # Cambiar según nombre de variable independiente
                   DFerror = df.residual(modelo.dcl), 
                   MSerror = dvmisc::get_mse(modelo.dcl),
                   alpha = 0.05,
                   group=TRUE,
                   main = NULL,
                   console=FALSE)
  return(comp[[5]] %>%
           rename("rdto" = "x$rdto") # Cambiar según nombre de variable respuesta
         )
}

Comparaciones de los niveles de B dentro del nivel A1

fcomp1(FA_A1)
       rdto groups
B1 5.334000      a
B2 4.885937      b

Comparaciones de los niveles de B dentro del nivel A2

fcomp1(FA_A2)
       rdto groups
B2 6.629250      a
B1 6.315125      a

Comparaciones de medias de los niveles de B dentro de cada nivel de A

FB_B1 <- data %>% filter(B=="B1")
FB_B2 <- data %>% filter(B=="B2")

fcomp2 <- function(x){
  comp <- LSD.test(x$rdto, # Cambiar según nombre de variable respuesta
                   x$A, # Cambiar según nombre de variable independiente
                   DFerror = df.residual(modelo.dcl), 
                   MSerror = dvmisc::get_mse(modelo.dcl),
                   alpha = 0.05,
                   group=TRUE,
                   main = NULL,
                   console=FALSE)
  return(comp[[5]] %>%
           rename("rdto" = "x$rdto") # Cambiar según nombre de variable respuesta
         )
}

Comparaciones de los niveles de A dentro del nivel B1

fcomp2(FB_B1)
       rdto groups
A2 6.315125      a
A1 5.334000      b

Comparaciones de los niveles de A dentro del nivel B2

fcomp2(FB_B2)
       rdto groups
A2 6.629250      a
A1 4.885937      b

Métodos Estadísticos para Experimentos Agrícolas en R